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Abstract 


In  this  research  project  modeling  methodologies  were  developed,  implemented  and  tested 
both  for  rapid  thermal  finite  element  analysis  of  small-scale  integrated  circuit  features  in  MCMS 
and  for  thermal  stress  finite  element  analysis  of  level  1  (chip-to-substrate)  wirebond  and  TAB  in¬ 
terconnects.  Due  to  the  small  size  of  IC  "microfeatures",  a  three  step  sequential  analysis  method¬ 
ology  was  developed  which  is  initiated  with  a  macroscopic  thermal  analysis  of  the  entire  MCM 
package.  The  macroscopic  finite  element  thermal  analysis  is  then  followed  by  two  successive  fi¬ 
nite  element  thermal  submodels  of  first  the  hottest  die  and  then  of  the  hottest  die  microfeature.  In 
this  manner  the  thermal  analysis  process  mathematically  "zooms"  into  the  hottest  IC  microfeature 
without  resorting  to  supercomputer-size  finite  element  models  of  the  MCM. 

In  a  similar  manner  a  two  step  sequential  thermal-stress  finite  element  submodeling  analysis 
procedure  was  developed  for  thermally  induced  stress  analysis  of  the  most  highly  stressed  wire- 
bond  or  TAB  interconnect  in  a  MCM  package.  For  automation  purposes  both  the  IC  thermal  sub¬ 
modeling  and  the  interconnect  elastostatic  submodeling  methodologies  were  implemented  into  an 
existing  blackboard-based,  object-oriented  MCM  software  design  tool  called  the  Intelligent  MCM 
Analysis  (IMCMA). 

Finally,  finite  element  algorithms  for  thermal  stress  analysis  involving  brick  and  tetrahedral 
elements,  as  well  as  state-of-the  art  finite  element  error  estimation  algorithms,  were  implemented 
into  FEECAP,  the  existing  finite  element  analysis  code  employed  by  IMCMA.  The  finite  element 
algorithms  were  validated  by  benchmark  comparisons  with  a  commercial  finite  element  code  and 
by  mesh  convergence  studies. 
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1  Introduction 


Multichip  modules  (MCMs)  are  high  performance  microelectronic  devices,  consisting  of 
several  chips  mounted  and  interconnected  to  a  multilayer  substrate  (Figure  1  [32]).  MCMs  are 
currently  used  in  military,  aerospace  applications  and  in  mainframe  computers  [1].  Figure  2  illus¬ 
trates  the  key  features,  including  interconnect  classification,  of  a  basic  low  cost  module  containing 
three  dies,  which  are  often  referred  to  as  chips.  On  the  top  surface  of  each  die  is  an  integrated  cir¬ 
cuit.  The  dies  may  be  mechanically  mounted  to  the  common  circuit  base  (substrate)  with  a  die  at¬ 
tach  adhesive  and  separate  first-level  electrical  interconnections  as  shown  in  the  figure,  or  in  the 
case  of  the  flip  chip  technology  a  solder  bump  array  is  used  to  both  mechanically  mount  the  dies 
and  function  as  first-level  electrical  interconnections.  Similarly,  the  MCM  package  itself  may  be 
mechanically  mounted  to  a  printed  wiring  board  by  an  adhesive  with  separate  second-level  inter¬ 
connects  for  electrical  connection,  or  the  package  may  be  both  electrically  and  mechanically  to  the 
printed  wiring  board  using  a  ball  grid  array  of  solder  bumps. 

Figure  3  shows  two  types  of  common  first  level  interconnects:  tape  automated  bonds  (TAB) 
and  wirebonds  [32].  MCMs  can  fail  through  the  fracture  or  debonding  of  a  first  level  interconnect 
such  as  the  wirebond  or  TAB  (tape  automated  bonding)  bond.  Wirebonding  is  the  most  common 
interconnect  technology  used  to  make  electrical  connections  between  the  chips  and  the  substrate. 
A  wirebond  is  a  wire  bonded  at  its  ends  to  the  substrate  and  a  chip,  by  ultrasonic,  thermocompres¬ 
sion,  or  thermosonic  bonding.  The  wire  material  is  usually  either  gold  or  an  aluminum- 
magnesium  alloy.  Tape  automated  bonding  can  produce  a  much  higher  interconnect  density  than 
wirebonds.  TAB  bonds  consist  of  patterned  metal,  usually  copper,  attached  to  polymer  tape  [2]. 
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Figure  1:  Low  Cost  MCM  (source:  Multichip  Module  Technologies  and  Alternatives:  The 
Basics,  edited  by  D.A.  Doane  and  P.D.  Franzon,  Van  Nostrand  Reinhold,  NY,  1993,  p.  80) 

Another  form  of  mechanical  failure  in  MCMs  is  thermal  derating-  the  overheating  of  a  com¬ 
ponent.  Both  temperature  and  stress  values  must  be  examined  by  the  MCM  designer  to  determine 
an  effective  mechanical  design. 
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PRINTED  WIRING  BOARD 

Figure  2:  MCM  Architecture  (schematic),  (adopted  from  Multichip  Module  Technologies  and 
Alternatives:  The  Basics,  edited  by  D.A.  Doane  and  P.D.  Franzon,  Van  Nostrand  Reinhold, 
New  York,  1993,  p  5.) 
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WIREBOND 
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Figure  3:  Common  types  of  first  level  connections:  chip  to  common  circuit  base  (adopted 
from  Multichip  Module  Technologies  and  Alternatives:  The  Basics,  edited  by  D.A.  Doane  and 
P.D.  Franzon,  Van  Nostrand  Reinhold,  New  York,  1993,  p  7.) 
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MCMs  are  extremely  expensive  and  time  consuming  to  manufacture,  so  it  is  important  for 
design  engineers  to  have  tools  for  determining  promising  designs  early  in  the  design  process.  Re¬ 
search  has  shown  that  evaluation  of  several  design  candidates  early  in  the  design  cycle  can  reduce 
the  cost  and  length  of  the  design  process  and  improve  the  robustness  of  the  MCM  [3].  Thus,  in 
the  domain  of  these  advanced  microelectronic  packages,  rapid  design  evaluation  via  modeling  and 
simulation  tools,  in  lieu  of  expensive  and  time-consuming  prototyping  and  testing  of  various  alter¬ 
native  package  designs,  is  critical  for  selecting  the  most  promising  design  for  subsequent  proto¬ 
typing  and  reliability  testing.  Note  that  given  the  current  limitations  in  predicting  MCM 
reliability,  the  purpose  of  a  present-day,  design-oriented  simulation  tool  should  not  be  to  eliminate 
MCM  prototyping  and  reliability  testing  but  rather  to  minimize  the  number  of  MCMs  prototyped, 
tested,  and  re-engineered.  It  should  also  be  emphasized  that  in  a  production  environment,  espe¬ 
cially  in  the  computer  industry,  time  to  market  is  absolutely  critical.  Therefore,  a  design-oriented 
MCM  simulation  tool  must  be  as  automated  as  possible  and  provide  feedback  in  minutes  or  no 
more  than  a  few  hours  of  real  time. 

The  Intelligent  Multicliip  Module  Analyzer  (IMCMA)  is  an  automated  finite  element 
analysis  (FEA)  based  mechanical  design  system  for  MCMs.  IMCMA  automatically  creates  a 
minimum  degree-of-freedom  finite  element  mesh  of  a  MCM  from  user-defined  high-level  infor¬ 
mation  and  computes  the  temperature  distribution  in  the  large-size  MCM  features  (such  as  sub¬ 
strate,  dies,  die  attach,  package,  etc.)  through  finite  element  analysis.  Thus,  a  design  engineer  can 
use  IMCMA  to  quickly  simulate  the  macroscopic  thermal  and  mechanical  behavior  of  various 
MCM  package  design  configurations.  While  IMCMA  provides  the  design  engineer  with  impor¬ 
tant  information  on  the  macroscopic  thermal  mechanical  behavior  of  a  MCM  design,  it  does  not 
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provide  any  data  on  the  thermal  or  mechanical  behavior  of  microscopic  MCM  features,  such  as 
power  dissipating  integrated  circuit  (IC)  features  and  first-level  interconnects  whose  behavior  are 
often  critical  to  the  reliability  of  a  MCM  package. 

2  Objectives 

The  objectives  of  this  research  project  are  to  develop,  implement,  and  test  methodologies' 
for  predicting  both  the  thermal  behavior  of  power-dissipating  integrated  circuit  features,  such  as 
field  effect  transitions  (FET),  and  the  mechanical  behavior  of  MCM  first-level  interconnects,  spe¬ 
cifically  wirebonds  and  TABs.  It  is  especially  importance  that  these  methodologies  be  computa¬ 
tionally  efficient  and  rapid  with  sufficient  accuracy  needed  to  compare  alternative  packaging 
design  concepts  in  terms  of  relative  mechanical  MCM  device  reliability.  The  methodologies  will 
be  instantiated  in  proof-of-concept  software  and  integrated  into  IMCMA,  a  prototype  MCM  pack¬ 
aging  design  system.  IMCMA  is  discussed  in  more  detail  in  Section  4. 

3  Proposed  Methodology 

The  finite  element  method  is  a  well  established,  proven  algorithmic  method  for  numerically 
solving  the  differential  equations  of  equilibrium  which  govern  the  behavior  of  physical  systems. 
Flowever,  there  are  two  main  obstacles  to  overcome  with  regard  to  the  application  of  the  finite 
element  method  as  a  design  tool  to  predict  the  thermal  and  mechanical  behavior  of  microscopic 
MCM  features  such  as  IC  die  features  and  first-level  interconnects.  First,  development  of  an 
MCM  finite  element  model  with  a  generic  finite  element  code  requires  significant  human  time,  fi¬ 
nite  element  modeling  expertise,  and  MCM  modeling  expertise.  Secondly,  the  sheer  number,  min¬ 
ute  sizes,  and  the  interaction  of  these  microscopic  features  with  the  macroscopic  MCM  behavior 
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precludes  a  "brute  force"  modeling  approach  in  which  all  macroscopic  and  microscopic  MCM 
features  are  geometrically  represented  in  a  single,  large  finite  element  model.  Such  a  brute  force 
model  would  consist  of  tens  of  thousands  of  finite  elements  and  nodes,  requiring  enormous  com¬ 
putational  resources  and  time  to  develop  and  analyze  and  thereby  obviating  the  model's  usefulness 
as  a  MCM  design  tool. 

Clearly,  a  much  more  intelligent  modeling  methodology  is  needed  which  can  effectively  re¬ 
solve  these  obstaeles.  Such  an  approach  lies  in  the  merging  of  three  key,  and  previously  isolated, 
technologies:  1)  high-level,  object-oriented  data  representation,  2)  the  blackboard-based  problem 
solving  paradigm,  and  3)  finite  element  submodeling  techniques.  The  first  two  technologies  have 
been  effectively  realized  in  IMCMA  for  macroscopic  thermal  finite  element  analysis  of  MCMs. 
IMCMA  takes  advantage  of  the  characteristics  of  the  MCM  domain  to  make  modeling  simplifica¬ 
tions  that  greatly  decrease  the  modeling  and  analysis  time,  while  still  providing  results  that  are 
useful  to  the  designer. 

Finite  element  submodeling  involves  a  multistep  finite  element  modeling  and  analysis  proc¬ 
ess.  In  the  initial  step  a  macroscopic  finite  element  analysis  of  the  entire  system  is  conducted.  The 
system  analysis  results  are  used  to  identify  critical  regions  or  features  of  the  model  which  must 
now  be  analyzed  in  more  detail  in  a  second  finite  element  model  of  only  the  critical  region.  The 
inherent  coupling  between  this  critical  region  and  the  initial  model  is  handled  by  imposing  bound¬ 
ary  conditions  on  the  critical  region  which  have  been  appropriately  obtained  from  the  initial 
analysis.  Essentially,  the  technique  represents  a  mathematical  "zooming"  into  a  critical  feature  of 
the  model. 
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An  important  characteristic  of  finite  element  submodeling  is  that  the  representation  of  the 
critical  feature  in  the  initial  finite  element  model  can  be  greatly  simplified.  Indeed,  if  the  critical 
feature  has  negligible  effect  on  the  thermal  or  mechanical  behavior  of  the  system  model,  then  the 
critical  feature  can  even  be  omitted  in  the  initial  system  model  but  then  modeled  in  detail  in  the 
subsequent  critical-feature  model.  For  example,  in  many  MCM  devices,  wirebonds  and  TABs 
provide  relatively  insignificant  heat  transfer  and  structural  stiffness.  Thus,  wirebonds  or  TABs  can 
be  neglected  in  an  initial  thermomechanical  finite  element  model  of  the  entire  MCM  system.  The 
system  results  are  then  used  to  identify  the  critical  interconnect  location  in  the  MCM  system. 
Boundary  conditions  at  this  location  are  extracted  from  the  system  analysis  results  and  then  im¬ 
posed  on  a  detailed  finite  element  model  of  the  interconnect  to  obtain  its  mechanical  behavior. 

Submodeling  thus  involves  modeling  simplifications  based  on  domain-specific  knowledge, 
i.e.  an  understanding  of  the  MCM  domain  and  the  types  of  simplifying  modeling  assumptions 
which  can  be  made.  The  system  architecture,  data  representation  scheme,  and  problem  solving  ap¬ 
proach  employed  in  IMCMA  is  well  suited  for  encapsulation  of  domain-specific  submodeling 
simplication  knowledge  and  automation  of  the  entire  submodeling  process.  In  section  4  we  first 
present  a  brief  description  of  the  IMCMA  software  system,  followed  by  details  on  thermal  and 
first-level  wirebond  and  TAB  interconnect  submodeling.  In  Section  5  benchmark  results  are  pre¬ 
sented  for  both  thermal  and  interconnect  submodeling.  Verification  of  the  thermal  stress  finite  ele¬ 
ment  analysis  algorithms  implemented  in  the  finite  element  solver  FEECAP  as  part  of  this 
research  project  is  contained  in  Section  5.  Conclusions  are  presented  in  Section  6.  Appendices  are 
attached  which  provide  details  on  new  error  analysis  algorithims  implemented  in  EEECAP  for 
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assessing  and  controlling  the  discretization  error  for  thermal  stress  problems,  as  well  as  opera¬ 
tional  details  related  to  thermal  and  interconnect  submodeling. 

4  The  Intelligent  Multichip  Module  Analyzer  (IMCMA) 

IMCMA  is  a  sophisticated  design  tool  for  rapid  thermal  finite  element  analysis  of  MCM 
package  designs.  IMCMA  employs  is  built  upon  a  blackboard  system  architecture  tightly  inte- ' 
grated  an  object-oriented  database  []-[].  The  object  oriented  database  and  blackboard  system  ar¬ 
chitecture  greatly  facilitates  seamless  integration  of  submodeling  software  that  was  developed 
under  this  project.  Thus,  IMCMA  was  chosen  as  the  software  vehicle  for  implementing  the  MCM 
submodeling  methodologies. 

4.1  The  Object-oriented  IMCMA  Database 

Conventional  relational  database  management  systems  (DBMSs)  are  not  well  suited  for  ap¬ 
plications  such  as  computer-aided  design.  Relational  DBMSs  are  not  capable  of  supporting  the 
complex  relationships  between  mechanical  parts.  In  a  relational  DBMS,  geometric  models  are 
stored  in  a  large  number  of  tables.  Explicit  links  between  tables  needed  to  represent  interrelations 
between  parts  do  not  exist.  The  manipulation  of  these  tables  then  requires  large  computer  pro¬ 
grams  which  have  long  execution  times  [4]. 

Extensive  research  has  been  done  recently  in  developing  object-oriented  databases  which 
are  capable  of  representing  complex  relationships.  IMCMA  uses  a  database  that  supports  relation¬ 
ships  between  physical  components,  and  also  between  components  of  the  finite  element  model 
and  mesh.  We  chose  an  object-oriented  database  to  fulfill  this  requirement.  Storing  physical  and 
finite  element  data  in  objects  provides  quick  access  to  data  describing  any  entity.  This  is  due  to 
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the  data  structure  itself,  and  also  to  the  fact  that  the  data  is  stored  in  memory,  not  on  disk.  The  da¬ 
tabase  groups  related  information  together,  and  this  reduces  data  access  time.  If  IMCMA  re¬ 
trieves  nodal  coordinates  and  nodal  constraints,  it  does  not  have  to  search  through  different  data 
files  to  compile  this  information.  Everything  is  contained  in  the  nodal  object.  The  elimination  of 
disk  searches  also  greatly  improves  data  access  time,  so  IMCMA  is  faster  than  a  relational 
DBMS-based  system  would  be. 

In  the  IMCMA  database,  each  physical  and  finite  element  entity  is  an  object.  An  object 
contains  information  that  defines  the  entity,  and  also  information  that  defines  the  relationship  of 
that  entity  to  others.  For  example,  a  nodal  object  would  contain  not  only  the  nodal  coordinates, 
loads,  and  constraints,  but  also  its  connectivity  to  other  nodes,  the  elements  to  which  the  node  be¬ 
longs,  and,  if  applicable,  the  geometric  surface  which  contains  the  node. 

The  objects  are  organized  into  a  hierarchy  of  abstractions  or  object  classes.  Figure  4  shows 
a  simplified  hierarchy  for  a  three  component  MCM.  The  different  object  classes  give  the  system 
an  ability  to  focus  on  specific  sets  of  data.  In  each  step  of  the  analysis,  IMCMA  can  move  up  and 
down  the  object  class  hierarchy  to  an  appropriate  level  of  abstraction.  Unnecessary  details  per¬ 
taining  to  other  object  classes  are  hidden  from  the  current  object  class,  a  property  called  informa¬ 
tion  hiding.  In  this  way  the  data  is  more  manageable,  and  the  analysis  is  faster  and  more  efficient 
[5].  The  combination  of  abstraction  and  a  blackboard  system  is  particularly  effective,  because 
blackboards  are  adept  at  moving  among  multiple  levels  of  abstraction  during  problem  solving  [6]. 

The  object-oriented  database  also  allows  for  an  excellent  user  interface.  IMCMA  produces 
plots  of  the  geometry  and  the  finite  element  meshes  on  the  screen  while  running  the  analysis.  The 
user  can  choose  entities,  such  as  lines,  surfaces,  nodes,  and  elements,  with  the  mouse,  and  all  of 
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Figure  4:  Example  Object  Hierarchy 


the  data  associated  with  the  entity  is  displayed. 


4.2  Blackboard  Systems  and  Knowledge  Sources 

IMCMA's  object-oriented  database  exists  in  the  environment  of  a  blackboard  system.  The 
blackboard  approach  is  a  flexible  artificial  intelligence  (AI)  problem-solving  technique  [7].  We 
used  the  framework,  a  toolkit  for  the  construction  of  high-performance  blackboard  appli-  , 

cations,  to  build  IMCMA.  GBB  supplies  the  necessary  object-oriented  blackboard  database. 

A  blackboard  system  uses  three  basic  components  to  solve  problems  [8]: 

(1)  A  global  database  (blackboard)  containing  input  data,  partial  solutions,  and  other  data  that 
are  in  various  problem-solving  states. 

(2)  Knowledge  sources  (KSs)  which  are  independent  modules  that  contain  the  knowledge 
needed  to  solve  the  problem,  and  that  can  be  widely  diverse  in  representation  and  in  infer¬ 
ence  techniques. 

(3)  A  control  mechanism,  that  is  separate  from  the  individual  KSs  and  that  makes  dynamic 
decisions  about  which  KS  is  to  be  executed  next. 

IMCMA's  knowledge  sources  carry  out  all  modeling  and  analysis  functions,  including  the 

definition  of  geometry,  meshing  of  the  model,  application  of  loads  and  constraints  to  the  finite 
element  model,  and  finite  element  analysis  of  the  model. 

I 

The  current  knowledge  sources  contained  in  IMCMA  are: 

•  input-model-ks:  defines  the  device,  components,  and  materials 

•  adjust-model-ks:  adjusts  chip  geometry  for  better  mapmeshes 

•  complete-model-ks:  creates  2D  and  3D  lines,  points,  and  surfaces 

•  generate-mm-regions-ks:  generates  mapmesh  regions 

•  find-symmetry-ks:  calls  a  CLIPS  program  for  model  simplification  based  on  symmetry 

•  generate-2d-mesh-ks:  creates  a  2D  meshes  using  FASTQ,  a  two-dimensional  mesh 
generation  tool  from  Sandia  National  Laboratory 

I  is  a  product  of  Blackboard  Technology  Group,  Amherst,  M A. 
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•  extrude-component-ks;  creates  3D  meshes  of  components  by  extrusion  of  2D  meshes, 
using  GEN3D,  a  mesh  extruder  from  Sandia  National  Laboratory 

•  combine-3d-meshes-ks:  combines  3D  meshes  using  GJOIN,  a  mesh  joiner  from  Sandia 

•  analyze-3d-mesh-ks:  performs  finite  element  analysis  of  model  using  FEECAP 

•  model-submodel-transition-ks;  identifies  component  to  be  submodeled,  resets  variables 
and  parameters 

The  blackboard  contains  all  the  data  needed  for  the  operation  of  the  knowledge  sources. 
IMCMA  stores  the  data  in  GBB's  object-oriented  database,  in  order  to  optimize  data  retrieval  for 
the  finite  element  analysis. 

5  Modeling  Simplification 

One  of  the  major  objectives  of  IMCMA  is  to  quickly  provide  designers  with  information 
that  will  help  them  to  make  intelligent  design  decisions.  Modeling  simplifications  are  tools  used 
by  IMCMA  to  construct  macroscopic  finite  element  models  of  MCMs  that  still  provide  meaning¬ 
ful  results.  Finn  et  al.  [9]  break  down  mathematical  modeling  simplification  into  two  categories, 
geometric  and  phenomena  simplification. 

The  four  types  of  geometric  simplifications  are  dimensional  reduction,  geometric  symmetry, 
feature  removal,  and  domain  alteration: 

Dimensional  reduction:  It  may  be  possible  to  reduce  the  degree  of  spatial  or  temporal 
analysis.  For  example,  a  3-D  model  may  be  reduced  to  a  2-D  model,  or  a  transient  analy¬ 
sis  may  be  reduced  to  a  steady-state  analysis. 

(2)  Geometric  symmetry:  The  finite  element  model  can  often  be  reduced  by  taking  advantage 
of  the  geometric  symmetry.  The  application  of  symmetry  boundary  conditions  on  the  re¬ 
duced  model  provides  the  same  physical  representation  of  the  system. 
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(3)  Domain  alteration:  Certain  aspects  of  the  physical  domain  can  be  changed  in  order  to 


simplify  the  analysis. 

(4)  Feature  removal:  Features  or  components  which  do  not  greatly  effect  the  macroscopic 
behavior  of  the  physical  system  can  be  eliminated  from  the  finite  element  model.  Cor¬ 
rectly  done,  the  removal  of  features  can  produce  a  greatly  simplified  finite  element  model 
while  maintaining  solution  accuracy.  For  example,  the  feature  may  be  a  possible  failure 
site,  and  the  maximum  stress  in  the  feature  will  need  to  be  calculated. 

The  use  of  submodeling  provides  the  ability  to  remove  features  at  the  macroscopic  or  large- 
scale  level  while  retaining  the  ability  to  perform  a  detailed  analysis  on  a  microfeature  or  small- 
scale  component  of  interest.  For  example,  the  microfeature  may  be  a  possible  interconnect  failure 
site,  such  as  the  heel/bonding  pad  interface  of  a  wirebond,  and  therefore  the  maximum  stress  in 
the  wirebond  heel  needs  to  be  computed.  A  separate  finite  element  submodel  analysis  of  the  wire- 
bond  is  carried  out  with  the  boundary  conditions  for  the  submodel  nodes  interpolated  from  a  pre¬ 
viously  conducted  finite  element  analysis  of  the  large-scale  components  that  comprise  the  MCM 
package.  The  submodel  analysis  results  provide  a  stress  or  temperature  distribution  of  a  feature 
that  is  likely  to  be  critical.  By  identifying  the  critical  small-scale  region  from  the  macroscopic 
analysis  results  and  subsequently  modeling  and  analyzing  only  the  microfeature  in  this  region,  a 
huge  amount  of  computational  effort  is  saved.  A  finite  element  model  of  a  MCM  package  with 
all  large-scale  and  small-scale  features  included  at  the  onset  is  highly  impractical  and  may  be  un¬ 
necessary,  depending  on  the  domain,  since  it  may  be  possible  to  obtain  all  necessary  information 
through  efficient  submodeling  techniques. 
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The  second  category  of  mathematical  modeling  simplification  is  phenomena  simplification. 
Phenomena  may  be  removed  or  reduced.  Phenomenon  removal  is  the  omission  of  an  entire  phe¬ 
nomenon  from  an  analysis.  A  decision  to  complete  only  a  thermal  analysis  is  phenomenon  re¬ 
moval,  because  stress  is  completely  ignored.  Phenomenon  reduction  is  the  removal  of  a 
component  of  a  phenomenon.  Ignoring  radiation  in  a  heat  transfer  analysis  is  an  example  of  phe¬ 
nomenon  reduction. 

5. 1  The  Application  of  Modeling  Simplification  in  IMCMA 

Modeling  simplifications  employed  for  developing  macroscopic  MCM  models  are  drawn 
from  knowledge  of  the  MCM  domain.  At  this  point,  it  is  important  to  understand  the  stage  of  the 
design  process  in  which  IMCMA  will  operate.  One  of  IMCMA's  purposes  is  to  provide  mechani¬ 
cal  analyses  of  MCMs  whose  electrical  layout  has  been  completed.  The  initial  design  that  IM¬ 
CMA  will  analyze  has  a  layout  which  includes  chip  and  interconnect  placement.  The  die  to 
substrate  interconnects  are  neglected  in  the  intial  macroscopic  analysis.  This  modeling  simplica- 
tion  is  based  on  the  fact  that  while  interconnect  materials  generally  have  high  thermal  conductiv¬ 
ity,  they  have  very  small  cross-sectional  area  (on  the  order  of  1000  pm^)  which  prevents 
interconnects  from  being  an  important  heat  path.  Interconnects  are  also  flexible,  so  while  the  dis¬ 
placement  of  a  substrate  or  chip  may  stress  an  interconnect,  the  displacement  of  the  substrate  or 
chip  is  not  greatly  affected  by  the  presence  of  interconnects.  Thus,  the  removal  of  the  intercon¬ 
nects  from  the  macroscopic  model  of  the  MCM  will  not  significantly  affect  the  macroscopic  be¬ 
havior  of  either  temperature  or  displacement.  This  feature  removal  provides  a  large  reduction  in 
computation  time.  Interconnects  are  failure  sites,  however,  so  they  cannot  be  ignored  altogether. 
IMCMA' s  submodeling  capability  provides  the  opportunity  to  model  critical  interconnects  later  in 
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the  analysis  process,  without  inclusion  of  interconnects  in  the  macroscopic  model.  Submodeling 
in  IMCMA,  which  will  be  discussed  in  the  next  section,  makes  the  removal  of  these  interconnect 
features  possible. 

The  macroscopic  finite  element  model  can  be  altered  to  include  the  effect  of  a  feature  called 
the  die  attach  without  geometrically  representing  and  meshing  the  die  attach  feature.  The  die  at¬ 
tach  is  the  layer  of  adhesive  between  a  chip  and  the  substrate.  The  layer  is  very  thin  compared  to 
the  thickness  of  the  chip.  This  small  die  attach  thickness  value  would  result  in  high  aspect  ratio 
finite  elements  that  model  the  die  attach  in  the  macroscopic.  High  aspect  ratios  can  introduce 
large  numerical  errors  into  the  finite  element  solution  due  to  the  finite  precision  with  which  digital 
computers  can  represent  numbers.  To  decrease  the  aspect  ratios  of  finite  elements  modeling  the 
die  attach,  a  much  finer  discretization  in  the  xy  space  of  the  die  attach  would  be  required.  This  re¬ 
sults  in  a  significant  increase  in  the  total  number  of  finite  elements  required  to  model  the  entire 
MCM  package  and  defeats  the  purpose  of  developing  a  simplified  efficient  macroscopic  model  to 
achieve  fast,  fairly  accurate  results. 

A  method  which  is  a  form  of  domain  alteration  has  been  developed  for  this  contract  to  in¬ 
clude  the  effect  of  the  die  attach  without  adding  the  die  attach  to  the  model.  In  this  approach  the 
thermal  conductivity  of  the  chip  is  changed  to  include  the  effect  of  the  die  attach  without  physi¬ 
cally  modeling  the  die  attach.  This  circumvents  the  problems  of  either  high  aspect  ratio  elements 
or  too  many  elements  for  efficient  macroscopic  analysis.  An  equivalent  thermal  conductivity  is 
calculated  by  the  following  formula: 


K,,  =  H- 


k\h 


k\h2+k2h\ 


(1) 
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where  /z,  is  the  die  attach  height,  hj  is  the  chip  height,  H  is  the  total  height  (/z,  +  h^).  A:,  is  the  die 
attach  thermal  conductivity  and  k2  is  the  chip  thermal  conductivity. 

The  critical  temperature  of  an  MCM  will  always  be  at  the  top  of  a  chip  at  the  location  of  the 
heat  producing  integrated  circuits.  By  using  an  equivalent  thermal  conductivity,  accurate  tem¬ 
perature  values  are  obtained  at  the  top  of  the  chip  with  the  use  of  a  macroscopic  model. 

IMCMA  also  employs  another  form  of  domain  alteration  called  xy-adjustment  is  used  to 
produce  simpler  mapmeshes.  If  sides  of  two  chips  are  nearly  collinear,  very  small  elements  will 
be  necessary  to  fill  in  the  space  between  the  two  sides.  The  macroscopic  solution  will  not  be  af¬ 
fected  significantly  by  a  slight  change  in  chip  size  or  location,  so  the  dimensions  of  the  chips  can 
be  altered  so  the  sides  are  collinear.  These  alterations  produce  a  much  simpler  mapmesh  with  a 
significant  reduction  in  total  number  of  finite  elements.  Therefore,  a  large  amount  of  computa¬ 
tional  effort  can  be  saved  with  almost  no  loss  in  accuracy.  The  engineer  has  control  over  the 
amount  of  xy-adjustment  in  chip  dimensions  that  can  be  introduced  into  the  solution;  IMCMA 
has  an  input  parameter  that  controls  how  much  the  dimensions  will  be  allowed  to  change.  IM¬ 
CMA  also  has  a  z-adjustment  capability  for  changing  chip  thicknesses  for  better  mapmeshes. 

IMCMA  utilizes  several  other  types  of  modeling  simplification.  Code  is  being  developed  to 
identify  geometric  symmetries  of  an  MCM,  and  simplifying  an  MCM  finite  element  model  ac¬ 
cording  to  these  symmetries  will  reduce  the  number  of  nodes  by  at  least  a  factor  of  two.  Domain 
alteration  is  also  used  through  the  assignment  of  a  general  set  of  material  properties  to  a  multi¬ 
level  substrate  which  contains  levels  of  wiring.  Other  features,  such  as  the  lid  which  covers  the 
MCM  package,  are  left  out  of  the  analyses,  because  they  are  not  critical  and  have  little  effect  on 
the  physical  behavior  of  the  other  components. 
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6  Implementation  of  Submodeling  in  IMCMA 

Two  types  of  submodeling  are  incorporated  into  IMCMA  .  thermal  IC  die  features  and  wire- 
bond  and  TAB  interconnects.  Thermal  die  feature  submodeling  is  used  to  find  the  on-die  tempera¬ 
ture  distribution  of  the  integrated  circuit  with  the  highest  surface  temperature.  Interconnect 
submodeling  is  used  to  identify  possible  interconneet  stress  failures  in  candidate  designs.  In  both 
cases,  submodeling  uses  the  results  of  the  IMCMA  maeroscopie  thermal  or  statie  analyses  to  get 
the  boundary  conditions  for  the  more  detailed  submodel. 

6.1  Thermal  Submodeling 

Thermal  submodeling  involves  several  separate  and  sequential  steps  as  illustrated  in  Figure 
5  below.  The  entire  process  is  fully  automated  in  IMCMA.  The  first  step  in  thermal  submodeling 
is  the  identification  of  the  hottest  chip.  This  is  done  by  examing  the  results  of  the  macroscopic 
MCM  thermal  analysis.  In  the  macroscopic  MCM  analysis  the  heat  producing  integrated  circuit 
features  of  each  die  are  approximated  by  a  uniform  surface  flux  applied  to  the  top  surface  of  the 
die.  An  average  power  dissipation  value  is  applied  to  each  surface  to  ensure  that  the  total  power 
generated  is  equal  to  the  actual  power  dissipated  on  the  chip.  The  use  of  flux  power  dissipation  al¬ 
lows  the  mesh  to  be  independent  of  the  geometry  of  the  heat  producing  regions. 

The  chip  with  the  highest  temperature  is  then  submodeled.  The  first  submodel  includes  the 
entire  chip  and  a  simplified  representation  of  the  heat  producing  regions.  The  user  can  also 
choose  to  have  the  die  attach  modeled  in  the  first  submodel.  The  chip-base  temperature  values 
obtained  from  the  macroscopic  MCM  thermal  analysis  are  applied  as  ehip  submodel  boundary 
eonditions.  Temperature  values  at  this  layer  are  interpolated,  because  the  submodel  has  a  higher 
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mesh  density  than  the  macroscopic  analysis.  The  interpolation  introduces  some  error,  but  the 
analyses  are  for  the  comparison  of  designs,  so  some  accuracy  can  be  sacrificed. 

The  second  submodel  contains  more  exact  heat  source  information  and  only  part  of  the  chip. 
The  amount  of  the  chip  that  is  modeled  is  dependent  on  the  mesh  density  of  the  first  submodel, 
but  the  maximum  thickness  is  one  quarter  of  the  total  chip  thickness,  and  the  maximum  top  area  is 
15  percent  of  the  total  chip  area.  Temperatures  are  interpolated  not  only  to  the  bottom  submodel 


nodes  in  this  case  but  also  to  all  nodes  on  the  sides  of  the  second  submodel.  The  results  of  the 
analysis  give  the  designer  an  estimate  of  maximum  package  temperature.  The  probability  of  ther¬ 
mal  derating  can  be  assessed  from  these  results,  and  a  better  thermal  design  can  be  obtained  by  se¬ 
lecting  the  design  with  the  lower  maximum  temperature. 

The  whole  thermal  submodeling  process  is  documented  on  the  workstation  screen  for  the 
designer  through  graphics  and  textual  descriptions.  Geometry,  finite  element  meshes,  and  thermal 
analysis  results  are  displayed  graphically  for  the  model  and  submodels.  The  text  takes  the  user 
through  all  knowledge  source  activations  and  describes  what  each  knowledge  source  is  complet¬ 
ing.  The  whole  submodeling  process  is  built  into  the  knowledge  sources  described  in  Section  4.2 
Blackboard  Systems  and  Knowledge  Sources.  The  changes  that  have  been  made  to  IMCMA  in 
order  to  implement  thermal  submodeling  are  detailed  in  Appendix  D:  Implementation  of  Ther¬ 
mal  Submodeling. 

6.1.1  Verification  of  Thermal  Submodeling 

A  verification  of  IMCMA's  thermal  submodeling  capability  was  completed  by  comparing 
IMCMA' ?,  finite  element  solutions  and  interpolation  values  to  an  ANSYS®^  analysis  of  the  same 
MCM.  The  model  used  for  verification  is  a  simple  two  chip  model.  Each  chip  has  nine  heat 
sources  mounted  on  the  top  surface.  Please  note  that  the  power  dissipation  rates  are  not  realistic 
for  an  MCM,  and  that  a  model  of  an  actual  MCM  would  likely  contain  a  greater  quantity  of 
smaller  heat  sources.  The  model  is  only  used  to  validate  IMCMA's  thermal  submodeling  proeess. 

The  model  consists  of  two  4x4x0.25  mm  silicon  chips  mounted  on  a  14x8x1  mm  aluminum 
substrate  (Figure  6).  A  prescribed  temperature  surface  of  20  “C  is  applied  to  the  bottom  surface  of 

^  ANSYS®  is  a  registered  commercial  product  of  Swanson  Analysis  Systems,  Inc.,  Houston,  PA. 
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Figure  6:  Model  for  Verification  of  Thermal  Submodeling 


the  substrate.  For  the  initial  macroscopic  analysis,  uniform  power  dissipations  are  applied  to  the 
chips.  The  power  dissipation  value  for  each  chip  is  equal  to  the  sum  of  the  power  dissipations  of 
the  chip  heat  sources. 

The  macroscopic  analyses  were  completed  on  ANSYS  and  IMCMA.  The  same  mesh  is 
used  for  each  analysis,  and  this  mesh  is  shown  in  Figure  7.  Figure  7  shows  IMCMA's  graphical 
user  interface,  which  displays  model  geometry  and  the  finite  element  mesh.  Both  the  ANSYS  and 
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Figure  7:  IMCMA  Graphics  Display  for  Macroscopic  MCM  Model 


IMCMA  analyses  produce  a  maximum  temperature  of  25.5  degrees  on  Chip  1.  Chip  1  is  identi¬ 
fied  as  the  hottest  chip  and  is  used  for  the  first  level  submodel. 

IMCMA  creates  the  first  submodel  geometry  and  applies  prescribed  temperature  values  to 
the  bottom  of  the  chip.  These  values  are  interpolated  from  the  macroscopic  model  results.  The 
ANSYS  and  IMCMA  interpolation  values  for  the  verification  model  were  found  to  agree  exactly 
to  two  decimal  places.  The  finite  element  mesh  and  geometry  for  the  first  submodel  are  shown  in 
Figure  8.  After  the  temperatures  are  prescribed,  IMCMA  examines  the  heat  sources  found  on  the 
chip  and  extrapolates  the  values  to  the  first  submodel  nodes  as  point  heat  sources.  ANSYS  does 
not  perform  this  operation,  but  IMCMA' s  point  heat  source  values  were  examined  and  found  to  be 
correct.  The  IMCMA  submodel  analysis  produces  a  maximum  temperature  of  1 12.8  °C  at  loca¬ 
tion  (4.00,  4.00,  1.25).  The  ANSYS  analysis  produces  a  maximum  temperature  of  1 13.3  °C  at  the 
same  point.  IMCMA's  result  differs  from  ANSYS  by  only  0.51  percent. 

The  second  submodel  contains  a  small  region  around  this  point  of  maximum  temperature. 
The  size  of  the  region  is  equal  to  the  size  of  9  first  submodel  elements  (3x3x1).  Prescribed  tem¬ 
perature  values  must  be  applied  to  the  front,  back,  left,  right,  and  bottom  faces.  The  values  are  in¬ 
terpolated  from  the  first  submodel  results.  The  major  difference  between  the  first  submodel  and 
the  second  submodel  is  the  representation  of  the  heat  sources.  In  the  first  submodel  the  heat 
sources  are  applied  to  the  nodes  as  point  heat  sources.  In  the  second  submodel  the  actual  geome¬ 
try  of  any  heat  source  within  the  critical  region  is  included.  The  physical  heat  source  is  meshed 
along  with  the  critical  part  of  the  chip,  and  a  heat  generation  rate  is  applied  to  the  volume  of  the 
heat  source. 
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Figure  8;  IMCMA  Graphics  Display  for  First  Submodel 
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The  ANSYS  model  used  for  verification  had  a  coarser  mesh  than  the  IMCMA  model,  but 
the  results  were  still  very  close.  Figure  9  shows  the  IMCMA  second  submodel  and  mesh.  The 
maximum  temperature  in  the  ANSYS  model  is  59.9  at  location  (2.67,  4.00,  1.30).  IMCMA's  re¬ 
sults  compare  well,  with  a  value  of  61.5  at  location  (2.75,  4.00,  1.25).  The  difference  in  tempera¬ 
ture  values  is  likely  the  result  of  a  slight  difference  in  location  and  the  difference  in  mesh  density. 

It  should  be  noted  that  the  first  submodel  analysis  has  a  maximum  temperature  nearly  twice 
as  high  as  the  maximum  temperature  of  the  second  submodel  analysis.  This  is  due  to  the  fact  that 
point  heat  sources  are  used  to  model  IC  power  dissipation  in  the  first  submodel,  while  the  physi¬ 
cal  IC  power  dissipating  component  is  modeled  in  the  second  submodel  as  a  three-dimensional 
solid  with  volumetric  heating.  High  temperatures  will  be  produced  at  the  points  where  the  heat 
sources  are  concentrated.  These  high  temperature  values  can  still  be  used  to  correctly  identify  the 
critical  region  of  the  submodel,  tlirough  the  additional  consideration  of  the  size  of  the  physical 
heat  sources.  The  maximum  nodal  temperature  on  an  element  which  contains  heat  sources  is  di¬ 
vided  by  the  area  of  the  heat  sources  contained  in  that  element.  The  element  which  has  the  high¬ 
est  temperature  to  area  ratio  is  identified  by  IMCMA  as  the  center  of  the  critical  region. 

6.2  Interconnect  Submodeling 

The  second  type  of  submodeling  methodology  developed  as  part  of  this  contract  is  first 
level  chip-to-substrate  connections.  The  different  types  of  first  level  connections  in  MCMS  are 
[cite  RL-TR  92-95]. 

•  wirebond 

•  tape  automated  bonding  (TAB) 

•  controlled  collapse  chip  connection  (C4) 

•  polymer  overlay/direct  metallization  on  die  (high  density  interconnection) 
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•  laser  pantography 

The  scope  of  this  contract  is  limited  to  wirebonds  and  TABs  which  are  established  first- 
level  interconnection  technologies.  The  C4  configuration  also  called  "flip-chip"  or  "solder-bump" 
bonding  technology  provides  the  highest  possible  number  of  I/O's  and  is  also  suitable  where  low 
inductance  is  an  important  consideration  in  the  design.  The  main  drawback  in  this  process  is  the 
requirement  for  a  full  area  array  of  first  level  connection  bonding  pads  on  the  chip  surface.  The 
vast  majority  of  chip  makers  design  and  manufacture  chips  with  bonding  pads  on  the  chip  surface 
peripheral.  IBM,  the  developer  of  the  C4  technology,  is  the  key  exception.  Polymer 
overlay/direct  metallization  was  developed  by  General  Electric  and  has  now  been  licensed  to 
Texas  Instruments.  Laser  pantography  was  developed  by  Lawrence  Livermore  National  Labora¬ 
tory.  This  technology  has  now  been  transferred  to  a  DoD  contractor. 

The  most  common  bond  type  used  is  the  wirebond.  Figure  10  shows  some  details  of  a  wire- 
bond  from  the  modeling  viewpoint.  Intensive  research  has  been  conducted  in  the  area  of  wire¬ 
bonding  to  determine  optimum  contours  of  the  wirebond  and  the  stress  distribution  in  them.  One 
of  the  most  important  reasons  for  the  popularity  of  this  type  of  bonding  is  the  low  cost  involved  in 
the  process.  The  automated  wirebonding  procedures  have  lead  to  fast  and  efficient  bonding  con¬ 
nections  as  compared  to  the  manual  bonding  that  used  to  be  previously  employed.  Some  advan¬ 
tages  of  wirebonding  are 

•  It  is  possible  to  use  various  metallizations. 

•  Adaptibility  to  wider  substrate  tolerances. 

•  Amenability  to  visual  inspection. 

•  Satisfactory  operation  characteristics  and  reliability. 
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Figure  10  :  Wirebond  Pro/Engineer  Solid  Modeling  Parameters 

The  TAB  technology,  however,  seems  to  be  taking  over  the  wirebonding  technology 
mainly  because  of  the  number  of  connections  that  can  be  made  to  a  chip.  Figure  1 1  shows  the  lay¬ 
out  of  TAB  tape  of  the  inner  lead  bond  and  outer  lead  bond.  The  inner  lead  bonding  consists  of  at¬ 
taching  the  unsupported  conductors  on  the  tape  to  the  die.  This  is  normally  done  using  the  process 
of  "gang  bonding"  i.e  simultaneous  bonding  of  the  leads  to  the  bond  pads.  The  usual  bonding 
methodology  used  is  the  thermocompression  (T/C)  bonding.  Thermocompression  bonding  is  dis¬ 
cussed  later  in  this  section.  Once  the  inner  lead  bonding  is  completed,  the  chip  with  the  leads  is 
placed  on  the  substrate  and  the  leads  are  thermosonically  or  ultrasonically  bonded  to  the  the 
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substrate.  These  two  bonding  techniques  are  also  discussed  in  brief  later.  Some  advantages  of  the 

TAB  over  wirebond  include  the  following  [cite  IPC-MC-790  and  Micro.  Pack.  HB,  p.  434]. 

•  improved  electrical  performance  at  high  frequencies 

•  ability  to  test  at  AC  speed  and  temperature,  and  bum  in  the  die  prior  to  final  assembly 

•  improved  heat  dissipation 

•  lower  bond  profiles  and  improved  component  reliability 

•  shorter  interconnect  lead  lengths,  resulting  in  higher  speed  performance 

This  report  discusses  the  modeling  and  analysis  techniques  of  the  most  commonly  used 
chip-level  bonding  technologies  -  wirebond  and  TAB. 
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6.2.1  Bonding  Methodologies  and  Bond  Types 

The  type  of  manufacturing  technique  to  be  employed  while  bonding  the  lead  to  the  chip  or 
substrate  affects  the  strength  of  the  bond  and  its  operational  characteristics.  Wirebonds  are  nor¬ 
mally  bonded  using  three  different  bonding  techniques  -  ultrasonic  (U/S)  bonding,  thermocom¬ 
pression  (T/C)  bonding  and  thermosonic  (T/S)  bonding.  In  ultrasonic  bonding,  the  wire  is 
clamped  onto  the  bond  pad  and  the  bonding  is  achieved  by  applying  ultrsonic  energy  normally  at 
ambient  temperature.  A  metallurgical  "cold  weld"  between  the  wire  and  the  bond  pad  results. 
Thermosonic  bonding  is  similar  to  U/S  bonding,  the  only  difference  being  the  slightly  elevated 
temperatures  of  the  component  used  in  the  process.  Ball  bonding  and  wedge  bonding  can  be  car¬ 


ried  out  using  U/S  and  T/S  methods.  Figure  12  shows  a  schematic  of  the  ball  and  wedge  bonds. 
Thermocompression  (T/C)  bonding  does  not  use  ultrasonic  energy.  The  bonding  is  achieved  by 
forcing  the  wire  on  the  bond  pad  and  using  thermal  energy  to  raise  the  temperature  (~  300°C)  at 
the  bond/surface  interface.  T/C  bonding  is  also  used  for  ball  and  wedge  bonding.  However  due  to 
the  long  bonding  times  and  the  high  temperatures  required  in  the  process,  it  is  becoming  less 
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popular  than  the  other  two  methods.  The  process  of  determining  the  critical  intercormect,  its  mod¬ 
eling  and  analysis  are  discussed  in  the  following  sectiion. 

6.2.2  Identification,  Modeling  and  Analysis  of  the  Critical  Interconnect 

Interconnect  submodeling  interfaces  the  commercial  solid  modeler/mesh  generator 
Pro/ENGINEER  with  IMCMA.  Interconnects  are  geometrically  complex,  unlike  rectilinear 
MCM  components,  and  cannot  be  modeled  and  meshed  automatically  by  IMCMA.  Parametric 
models  of  intercormects  are  created  in  Pro/ENGINEER,  and  IMCMA  provides  values  for  the  nec¬ 
essary  parameters,  fires  Pro/ENGINEER  to  mesh  the  model,  imports  the  mesh  data  back  into  the 
IMCMA  database,  and  performs  the  finite  element  analysis.  Finally  the  results  are  also  read  back 
into  IMCMA. 

The  first  step  in  interconnect  submodeling  is  the  identification  of  a  critical  interconnect. 
Research  data  has  shown  that  critical  wirebonds  and  TAB  bonds  are  found  at  the  comers  of  the 
chips.  Displacements  and  temperatures  are  found  at  these  comer  locations  from  the  macroscopic 
thermal/static  analysis  of  the  MCM.  Displacement  causes  stress  in  the  material  by  bending  the  in¬ 
terconnect,  and  temperature  causes  stress  through  the  mismatch  in  coefficient  of  thermal  expan¬ 
sion  between  the  interconnect  and  the  chip  or  substrate.  Relative  importance  of  temperature  and 
displacement  would  have  to  be  determined  through  an  extensive  study,  so  until  this  can  be  done,  a 
simple  methodology  is  used  to  identify  the  critical  interconnect.  The  interconnect  at  a  corner  lo¬ 
cation  with  the  highest  relative  displacement  between  its  ends  is  taken  to  be  critical.  IMCMA 
checks  the  relative  displacements  at  all  comer  locations,  and  finds  the  maximum. 

The  interconnect  model  is  interactively  created  by  the  user  in  Pro/Engineer  and  meshed. 
The  displacements  at  both  ends  of  the  interconnect  are  applied  as  boundary  conditions  to  the 
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Pro/engineer  finite  element  mesh.  IMCMA  performs  a  finite  element  analysis  of  the  interconnect, 
and  the  designer  has  access  to  all  results,  which  are  stored  in  the  database. 

6.2.2.1  Identification  of  the  Critical  Interconnect 
In  order  to  zoom  into  the  critical  regions  of  the  MCM  and  to  extend  the  analysis  to  a  micro¬ 
scopic  level,  it  is  necessary  to  identify  areas  of  the  MCM  that  are  prone  to  failure.  Typically,  the 
solution  to  this  is  submodeling.  Submodeling  involves  determining  the  features  that  may  be  criti¬ 
cal  to  the  design,  modeling  them  and  carrying  out  the  analysis.  Boundary  conditions  on  the  sub¬ 
model  are  specified  using  information  from  the  macroscopic  analysis. 

The  interconnect  is  one  such  component  of  the  MCM  which  tends  to  fail  and  hence  needs  to 
be  submodeled.  Each  MCM  may  contain  anywhere  between  100  to  500  interconnects  and  analyz¬ 
ing  all  of  them  would  take  up  a  large  amount  of  computation  time.  This  necessitates  the  identifica¬ 
tion  of  the  critical  interconnect  given  a  set  of  interconnects. 

In  order  to  identify  the  critical  interconnect  from  a  user  input  of  interconnects  a  critical  in¬ 
terconnect  knowledge  source  was  developed.  This  knowledge  source  identifies  the  critical  inter¬ 
connect  based  on  the  largest  relative  displacement  between  the  start  and  end  point  of  the 
interconnect  and  stores  the  attributes  of  the  critical  interconnect  on  the  blackboard.  The  basic 
concept  for  the  identification  of  the  critical  interconnect  is  discussed  below. 

Initially  a  macroscopic  mesh  analysis  is  performed  with  only  the  substrate  and  chip  configu¬ 
ration  i.e.  excluding  the  intereonnects.  The  nodal  displacement  and  nodal  temperature  information 
of  the  macroscopic  mesh  is  then  present  on  the  blackboard  and  can  be  retrieved  when  required. 
The  interconnect  geometry  is  defined  by  its  start  and  end  points.  These  points  are  located  with  re¬ 
spect  to  the  macroscopic  mesh  configuration.  The  elements  containing  the  start  and  end  points  of 
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the  interconnect  are  identified,  and  the  displacements  are  interpolated  from  the  nodes  of  these  ele¬ 
ments  to  the  start  and  end  points.  The  relative  displacement  is  then  computed  for  the  interconnect. 
This  procedure  is  carried  out  for  each  "chip  corner  interconnect"  and  the  interconnect  with  the 
largest  relative  displacement  is  determined.  This  interconnect  is  stored  on  the  blackboard  for  fur¬ 
ther  analysis  ("zooming"  onto  the  critical  interconnect).  All  other  interconncts  are  deleted  from 
the  blackboard  database. 

A  brief  description  of  each  module  or  knowledge  source  which  was  developed/modified  for 
the  identification  of  the  critieal  interconneet  is  contained  in  Appendix  F,  Determination  of  Critical 
Interconnect . 

6.2,2.2  Modeling  Of  the  Interconnect 

As  discussed  earlier,  interconnects  are  sensitive  eomponents  of  an  MCM.  Size  and  con¬ 
struction  are  important  factors  in  interconnect  failure,  so  an  authentic  model  and  an  accurate 
analysis  procedure  are  necessary.  The  complicated  construction  features  of  interconnects  intro¬ 
duce  modeling  complexities,  which  if  neglected  may  lead  to  errors  in  the  final  results. 

The  two  classes  of  interconnects  that  have  been  our  main  area  of  concentration  in  this  report 
are  wirebonds  and  TAB  (tape  automated  bond)  bonds.  While  the  latter  does  not  pose  too  many 
modeling  challenges,  the  contour  of  the  wirebond  demands  the  use  of  a  variety  of  modeling  skills. 
A  basic  assumption  of  a  skilled  modeler  is  made  while  designing  and  executing  the  system.  The 
modeling  of  the  interconnect  is  carried  out  using  the  parametric  modeler  in  Pro/ENGfNEER. 

CONSTRUCTIONAL  AND  FUNCTIONAL  DETATT.S  OF  THE  TNTEF CONNECT- 

The  interconnect  consists  of  a  gold  or  aluminum  alloy  lead,  typically  about  0.005  inches  in 
diameter  (for  a  wirebond)  or  width  (for  a  TAB).  The  entire  span  of  the  interconnect  is  of  uniform 
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cross-section  except  the  ends  are  flattened  out  or  given  a  ball  shape  to  facilitate  bonding.  The  flat 
end  region  is  known  as  the  thermocompression.  It  is  the  region  on  the  interconnect  which  estab¬ 
lishes  the  electrical  connection  between  the  chip  and  the  substrate  and  also  serves  as  a  reinforce¬ 
ment  to  the  chip.  The  interconnect  has  a  sudden  change  in  the  cross-sectional  area  in  the  region  of 
the  thermocompression.  This  area,  referred  to  as  the  bond  heel,  is  the  most  critical  region  of  the 
intercormect  due  to  a  concentration  of  stresses.  Some  previous  tests  to  determine  fatigue  life  of  in¬ 
terconnects  have  shown  failure  at  the  bent  arc  (the  topmost  region  of  the  interconnect)  but  this  has 
been  mainly  due  to  surface  asperities,  manufacturing  flaws  and  material  defects.  Failure  normally 
occurs  at  the  heel  of  the  interconnect. 

MOnF.T  .TNG  CONSIDERATIONS  AND  MESHING: 

The  modeling  methodology  consists  of  interactive  use  of  Pro/ENGINEER.  The  critical  in¬ 
terconnect  needs  to  be  modeled  and  analyzed.  The  user  is  initially  queryed  by  IMCMA  for  the 
material  of  the  interconnect.  This  material  and  its  properties  must  be  defined  in  the  example  file. 
Depending  on  the  number  of  chips  in  the  model,  the  four  comer  locations  of  each  chip  are  con¬ 
sidered  to  be  the  end  points  of  the  interconnect.  User  input  for  the  span  of  the  interconnect  is 
taken  and  the  start  points  are  determined.  The  assumption  made  here  is  that  the  material  and  span 
of  all  four  intercoimects  of  a  chip  are  the  same. 

Once  the  interconnect  is  identified,  a  Pro/ENGINEER  window  is  started  up.  The  start  point 
indicates  a  point  on  the  substrate  while  the  end  point  is  located  on  the  chip.  Figure  10  on  page  27 
shows  the  start  and  end  points  of  the  interconnect  and  the  thermocompression  region.  To  start  the 
modeling  procedure,  the  user  needs  to  read  in  a  file  of  points  {points.dat)  that  define  the  intercon¬ 
nect  contour.  This  file  is  created  by  IMCMA  while  determining  the  location  of  the  critical 
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interconnect.  It  consists  of  five  points  in  case  of  a  wirebond  (the  start-point,  end-point  and  three 
other  intermediate  points)  and  four  points,  the  start  and  end-points  and  two  other  intermediate 
points  in  case  of  a  TAB  bond.  These  points  are  read  in  as  datum  points  and  are  reference  locations 
for  the  interconnect  model.  The  start  point  is  always  at  the  coordinates  (1,1,1)  and  the  end  point  is 
displaced  relatively  from  the  start  point  depending  on  the  span  of  the  interconnect.  These  start  and 
end  point  coordinates  may  not  match  those  of  the  critical  interconnect.  Only  the  relative  x,  y  and  z 
offsets  are  the  same.  This  change  is  done  to  simplify  the  task  of  modeling  and  also  the  analysis 
procedures.  From  here  on  the  modeling  is  carried  out  by  a  skilled  modeler.  There  may  be  a  large 
number  of  methods  that  may  be  employed  in  modeling  the  interconnect.  The  choice  of  this 
method  is  left  to  the  user.  Some  methods  are  suggested  in  Appendix  G. 

Once  the  model  has  been  created,  the  user  meshes  the  model  using  the  tet-mesh  facility  in 
the  FEM  module  of  Pro/ENGINEER.  The  default  value  of  global  element  length  may  not  be  used 
as  the  resulting  meshes  have  a  large  percentage  of  distorted  elements.  Once  the  automatic  mesh 
has  been  created,  the  user  outputs  the  model  through  Pro/Engineer  to  an  ANSYS  input  file  (log 
file)  using  the  interfacing  capability  of  Pro/ENGINEER.  The  file  must  be  named  bond.ans.  The 
element  type  used  is  a  linear  4  noded  tetrahedral  element.  At  this  point  there  are  no  boundary  con¬ 
ditions  applied  to  the  model.  The  model  may  be  stored  by  the  user  and  he/she  may  exit 
Pro/ENGINEER. 

6.2.2. 3  Reading  Mesh  Data 

Once  the  ANSYS  input  file  of  the  interconnect  mesh  has  been  created,  the  mesh  data  is  read 
into  IMCMA  from  the  bond.ans  file.  All  nodes  and  elements  in  the  model  are  read  in  and  created 
as  objects  on  the  blackboard.  All  nodes  of  the  interconnect  are  created  as  instances  of  the  unit 
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class  intercoviTicct-Sd-nodc  while  the  elements  are  created  as  instances  of  the  unit  class 
interconnects d-element .  Each  instance  of  the  node  contains  information  on  the  x-y-z  coordinates 
and  is  linked  to  the  respective  elements  through  a  link  slot.  The  name  assigned  to  the  node  is  its 
number.  The  instances  of  the  elements  hold  information  on  the  nodal  instances  that  make  up  the 
element  and  the  nodal  connectivity  represented  as  a  node-list.  The  mesh  data  is  linked  together  to 
enhance  information  retrieval  using  different  GBB  commands. 

6.2.2.4  Analyzing  The  Model 

The  blackboard  contains  the  following  information  at  this  stage  of  the  modeling  process. 

•  Information  about  the  nodes  and  elements. 

•  Information  about  the  material  of  the  interconnect  and  its  properties  which  are  accessed 
from  the  material  library  . 

•  Start  and  end  displacements  of  the  critical  interconnect  obtained  from  the  previous  knowl¬ 
edge  sources. 

The  start  and  end  displacements  of  the  critical  interconnect  are  considered  to  be  uniform 
over  all  the  nodes  that  are  in  contact  with  the  substrate  and  the  chip  respectively.  These  are  the 
only  constraints  applied  to  the  model  during  the  analysis. 

The  nodes  that  lie  at  the  start  and  end  are  determined  using  their  z-coordinates  and  the 
length  (x-distance).  In  case  of  a  wirebond,  this  x-distance  is  the  thermocompression  length.  The 
nodes  on  the  interconnect  along  the  surface  of  contact  between  the  interconnect  and  the 
substrate/chip  are  assigned  the  start  displacements  or  end  displacements  depending  on  whether  the 
node  lies  at  the  start  or  end  of  the  interconnect. 
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The  mesh  information  is  now  written  into  the  INPUT.DAT  file.  This  file  is  the  input  for 
IMCMA's  finite  element  analysis  eode  FEECAP.  The  data  written  out  to  this  file  includes  nodal 
data,  elemental  data,  material  data,  loads  and  other  boundary  conditions  for  the  analysis.  Addi¬ 
tionally  it  contains  flags  for  the  analysis  type  and  the  type  of  element  being  used.  In  this  case  the 
analysis  is  always  a  static  analysis  (flag  lATYP  =  2)  and  the  element  is  always  a  4  noded  tetrahe¬ 
dral  element  (NELTYP  =  3).  The  analysis  is  carried  out  and  the  mesh  results,  namely  nodal  dis¬ 
placements  and  stresses,  are  calculated. 

The  analysis  results  are  written  out  to  the  file  RESULTS.OUT.  In  addition  the  displace¬ 
ments  and  stresses  are  written  to  temporary  files  DISP.OUT  and  STR.OUT  to  facilitate  reading 
the  results  into  IMCMA. 

Reading  the  analysis  results  constitutes  assigning  the  nodal  displacements  (x,  y,  z)  and 
stresses  (x,  y,  z,  xy,  yz,  zx,  and  Von-Mises  stress)  to  the  corresponding  accessors  of  the  slots  in 
the  unit  class  interconnect-3d-node.  The  displacements  and  stresses  are  stored  as  a  list  of  three  and 
seven  elements  respectively.  This  completes  the  entire  run  of  the  interconnect  submodeling  mod¬ 
ule  m  IMCMA. 

7  Thermal  Stress  Formulation  and  Benchmarking; 

7.1  The  Eight  Noded  3D  Linear  Brick  Element 

The  eight  noded  linear  brick  element  is  probably  the  most  widely  used  element  type  in  a  va¬ 
riety  of  design  engineering  problems  that  require  an  accurate  3D  finite  element  meshing  and 
analysis.  Within  IMCMA,  since  the  domain  is  restricted  to  MCM  package  design,  the  brick  ele¬ 
ment  has  been  chosen  for  macroscopic  finite  element  analysis  of  the  MCM.  The  formulation  and 
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system  equations  for  the  eight  noded  thermal  brick  element  were  included  in  the  June  1993  tech¬ 
nical  report  [31].  The  following  section  lists  some  of  the  element  formulation  equations  for  the 
eight  noded  brick  element  for  static  analysis.  A  more  detailed  description  of  a  generic  3-D  ele¬ 
ment  formulation  for  static  analysis  can  be  found  in  Appendix  E. 

7.1.1  Element  Formulation 

The  formulation  of  the  static  element  involves  the  applied  static  loads  and  initial  strains  due 
to  applied  thermal  loads.  The  following  are  the  element  equations  for  the  eight  noded  isoparamet¬ 
ric  brick  element.  The  equilibrium  equations  for  static  analysis  case  for  element  e  are  given  by: 


[K‘]{d‘}  =  {r\)  +  {r}]  +  {rl}  (2) 

where 

[K^]  -  (3) 

{r\}  =  ^WYmdT 

{r}]  =  \[N‘Y{f)da  (5) 

{n.}  =  \[B‘Y{E‘]Yl}dO.  (6) 


and  Q®and  T^"  donate  the  element  volume  and  element  surface  area,  respectively.  For  static 
analysis  the  matrix  [K*"]  is  the  element  stiffness  matrix,  {d*" }  is  the  vector  of  element  nodal  dis¬ 
placements,  [F]  is  the  material  constitutive  matrix  of  the  element,  [B"]  is  the  element  strain  dis¬ 
placement  matrix,  and  },  {r^},  and  {rl^ }  are  the  element  load  terms  due  to  applied  surface 
tractions  {f  },  applied  body  loads  {f},  and  initial  strains  {eq},  respectively.  The  elemental  ini¬ 
tial  strains  {Sq}  are  obtained  from  a  finite  element  thermal  analysis.  A  detailed  derivation  of  the 
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element  equations  for  the  eight  noded  isoparametric  elastostatic  brick  element  is  presented  in  Ap¬ 
pendix  E. 

7.1.2  Benchmarking  of  Results 

To  validate  the  results  obtained  by  FEECAP  an  example  problem  was  selected  and  the  ther¬ 
mal  and  stress  analysis  results  were  benchmarked  with  ANSYS  (version  5.0).  The  example  cho¬ 
sen  was  a  single  chip  surface  mounted  on  a  substrate  and  subjected  to  thermal  static  loads  and 
constraints.  Figure  13  illustrates  the  location  and  magnitude  of  the  applied  thermal  loads  and 
boundary  conditions.  The  element  types  chosen  in  ANSYS  for  thermal  analysis  and  static  analysis 
were  8-noded  brick  elements  of  type  SOLID70  and  SOLID45  respectively. 

Tables  1,2  and  3  show  the  nodal  values  obtained  from  ANSYS  and  the  percent  error  (%Eit) 
of  the  results  computed  by  FEECAP  with  respect  to  the  ANSYS  solution.  The  nodal  temperature 
values  were  found  to  be  accurate  upto  six  significant  digits  and  the  %Err  values  in  the  flux  com¬ 
ponents  were  less  than  0.04%.  The  %Err  values  in  the  displacement  components  were  within 
0.004%  and  for  the  cauchy  stress  components  the  maximum  %Err  value  was  within  0.005%. 
These  low  error  values  validate  the  finite  element  computations  performed  by  FEECAP  when 
compared  with  ANSYS. 
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Uniform  thickness  =  1.0 

Thermal  loads:  applied  flux,  ql  =  10  on  top  surface  of  element  10  (chip) 
applied  flux,  q2  =  -8  on  top  surface  of  element  1 
heat  generation  Q  =  5  by  element  7 
point  heat  sources  q  =  6  at  nodes  3  and  1 9 
prescribed  T  =  30  oif  right  face  of  substrate 
all  other  surfaces  are  insulated 

Static  constraints:  x,  y  and  z  constraints  on  bottom  surface  of  substrate 

Figure  13:  Example  of  a  chip  surface  mounted  on  a  substrate  for  benchmarking  the 
combined  thermal-static  stress  computations. 
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KEY :  T  -  temperature,  q^,  Qy  and  =  thermal  flux  components 

. ,  „  (ANSYS-value)  -  (FEECAP-value) 
(ANSYS-value) 


Table  1.  Thermal  Solution  of  the  Combined  Thermal  Static  Stress  Analysis  for  the 
Benchmark  Example 


40 


U-x  (xlO) 

%Err 
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KEY:  U^,  Uy  and  =  displacement  components 

(ANS YS-value)  -  (FEECAP-value)  ,  ^  „ 

%  Err  =  - - -  *  1  oO 

(ANSYS-value) 


Table  2.  Displacement  Solution  of  the  Combined  Thermal  Static  Stress  Analysis  for 
the  Benchmark  Example 
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KEY :  a^,  a^,  a^,  a^y,  and  =  Cauchy  stress  components 


%  Err  = 


(ANSYS-value)  -  (FEECAP-value) 
(ANSYS-value) 


*100 


Table  3.  Static  Stress  Solution  of  the  Combined  Thermal  Static  Stress  Analysis  for 
the  Benchmark  Example 
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7.2  The  Four  Noded  Linear  Tetrahedral  Element: 


Interconnects  form  a  class  of  objects  that  have  a  complex  shape  and  contour.  The  develop¬ 
ment  of  a  good  mesh  of  this  component  poses  a  problem.  A  poor  mesh  will  affect  the  results  ob¬ 
tained  from  the  finite  element  analysis.  The  four  noded  tetrahedral  element  was  found  to  be  one  of 
the  element  types  suitable  for  mesh  generation.  This  element  is  most  commonly  used  in  the  auto¬ 
matic  mesh  generators  of  commercial  software.  Additionally  this  element  type  was  one  of  the  few 
types  supported  by  the  finite  element  preprocessor  in  Pro/ENGINEER.  The  new  version  of  FEE- 
CAP  (ver.  2.5)  has  been  suitably  modified  to  provide  the  capability  to  carry  out  finite  element 
analysis  using  meshes  made  up  of  four  noded  linear  tetrahedral  elements.  Fifteen  new  subroutines 
have  been  added  and  changes  have  been  made  to  the  existing  FEECAP  routines  to  accomplish  this 
task. 

7.2.1  Element  Formulation: 

The  four  noded  linear  tetrahedral  element  formulation,  unlike  the  eight  noded  linear  brick 
element  formulation,  does  not  require  use  of  Gaussian  quadrature  to  integrate  element  matrices 
and  load  vectors.  Closed  form  expressions  for  these  quantities  are  obtained  without  employing  nu¬ 
merical  integration  or  mapping  of  elements  from  the  real  to  the  natural  coordinate  space.  All  cal¬ 
culations  are  done  at  nodal  points  of  the  element,  which  reduces  the  task  of  extrapolating  results 
from  the  gauss  points  as  in  case  of  an  isoparametric  element.  The  formulation  involves  the  consid¬ 
eration  of  different  cases  of  boundary  conditions  for  thermal,  static  and  combined  thermal  and 
static  analysis.  Figure  14  shows  the  nodal  numbering  scheme  followed  for  the  tetrahedral  element. 
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Figure  14:  The  Tetrahedral  Element 

The  following  is  a  brief  summary  of  the  formulation  of  element  equations  and  the  handling 
of  various  boundary  conditions  in  FEECAP.  Additional  mathematical  details  are  given  in  Appen¬ 
dix  I. 

7.2.1. 1  Thermal  Analysis 

The  mathematical  formulation  of  the  finite  element  method  involves  the  transformation  of 
governing  partial  differential  equations  and  boundary  conditions  of  a  physical  system  into  a  set  of 
element  algebraic  equations.  The  element  equations  are  assembled  into  the  system  equations  and 
then  solved  simultaneously. 

For  a  three  dimensional  thermal  analysis  the  following  element  equations  may  be  written 

[i^^]{r}  =  {4}  +  Kg}  +  {r^}  (7) 

where 
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(8) 

(4)= 

Ft’ 

(9) 

1 

{rg}  =  \{WVQ’^da 

(10) 

{rj,}  =  IvmVh'^ndr 

(11) 

In  the  above  equations  [K"]  is  the  element  conductance  matrix,  is  the  element  material 
conductivity  matrix,  [If]  is  the  matrix  of  shape  functions,  {F}  is  the  vector  of  element  nodal  tem¬ 
peratures,  and  the  {/}  vectors  are  the  element  nodal  heat  loads  due  to  different  loading  condi¬ 
tions.  The  first  load  term,  is  the  due  to  prescribed  flux  ^^on  element  face(s),  {rg}  is 

due  to  internal  volumetric  heat  generation  Q"  within  an  element,  and  is  due  to  a  convection 
boundary  condition  with  convection  coefficient  h  and  ambient  temperature  T^.  All  the  above 

equations  are  written  at  the  elemental  level  (hence  the  superscript  e). 

7.2.1.2  Static  Analysis 

The  governing  equations  for  the  static  analysis  with  the  tetrahedral  element  are  the  same  in¬ 
tegral  equations  developed  for  the  brick  element  listed  in  Equations  (2)  -  (66)  on  37page  37. 
FEECAP  also  has  the  capability  to  handle  prescribed  nodal  displacements  for  the  tetrahedral  ele¬ 
ment.  The  load  vector  is  adjusted  due  to  these  prescribed  displacements,  if  any,  before  computing 
the  other  nodal  displacements  in  the  model. 

The  formulation  uses  natural  coordinates  (area  and  volume  coordinates).  All  calculations 
and  evaluations  of  integrals  are  done  using  these  coordinates.  The  area  of  an  element  face  m 
terms  of  x  and  y  coordinates  of  its  nodes  is  given  by 
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1  1  1 

A'^  =  ^det  xi  X2  X3  (12) 

y\  yi 

and  the  volume  of  an  element  in  terms  of  the  coordinates  of  its  nodes  is  given  by 

1  x\  y\  z\ 

1  X2  yi  Z2 
1  X3  y3  Z3 
1  X4  y4  Z4 

Integration  of  a  polynomial  in  area  coordinates  over  the  triangle  area  is  accomplished  by  the 


following  formula 


k\l\m\ 

(2  +  A:  +  /  +  m)! 


(14) 


where  A  is  the  entire  area  of  the  triangular  face,  k,  1,  m  are  non-negative  integers  and  the  con¬ 
straint  given  by  the  following  equation  is  satisfied. 


4i+^2+^3  =  1 


Similarly  for  the  volume  coordinates 


A  more  detailed  treatment  of  related  equations  and  formulae  may  be  found  in  Appendix  I.  A 
generic  formulation  for  a  three-dimensional  element  is  given  in  Appendix  E. 


7.2.2  Benchmarking  of  Results: 

All  results  for  the  tetrahedral  elements  have  been  successfully  benchmarked  with  ANSYS 
5.0.  Three  examples  have  been  presented  in  this  report,  one  each  for  the  thermal,  static  and  com¬ 
bined  analysis.  The  model  used  for  all  the  examples  is  a  cube  of  size  2x2x2  inches  (see  Figure  15) 
with  thermal  conductivity  ^  =  10  Btu/hr  T  ,  coefficient  of  thermal  expansion  a  =  6.4  (10‘^)  units?, 
modulus  of  elasticity  E  =  54  (10^)  psi,  and  Poisson's  ratio  v  =  0.22.  Different  boundary  conditions 
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were  imposed  on  the  model  for  each  analysis  case.  The  results  were  compared  with  those  obtained 
from  ANSYS.  The  corresponding  element  type  used  in  ANSYS  was  the  collapsed  brick  element 
(SOLID70  for  thermal  analysis  and  SOLID45  for  static  analysis). 


Figure  15:  Thermal,  static,  and  combined  thermal-static  benchmark 
example  for  tetrahedral  element 


7.2.2.1  Thermal  Analysis  Results  for  Tetrahedral  Element 

The  following  are  the  comparison  results  for  thermal  analysis.  The  input  model  consists  of 
prescribed  temperatures  at  nodes  4,  6,  8  (refer  to  Figure  15  for  details)  and  flux  boundary  condi¬ 
tions  are  imposed  on  two  element  faces.  Table  4  shows  the  temperature  solution  for  the  above 
problem.  As  can  be  seen,  the  FEECAP  nodal  temperatures  have  five  significant  digit  accuracy 
compared  to  ANSYS  results.  Table  5  shows  the  benchmarked  flux  results. 


NODE 

TEMPERATURE 

ANSYS 

FEECAP 

1 

89.463 

89.463 

2 

92.301 

92.301 

3 

80.891 

80.891 

4 

40 

40 

5 

96.087 

6 

80 

7 

90.373 

90.373 

8 

90 

90 

9 

78.319 

78.319 

Table  4.  Temperature  Solution  Benchmark  Comparison  for  Thermal  Analysis 

(tetrahedral  element) 
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FLUX-X 

FLUX-Y 

FLUX-Z 

TOTAL  FLUX 

NODE 

ANSYS 

FEECAP 

ANSYS 

FEECAP 

ANSYS 

FEECAP 

ANSYS 

FEECAP 

1 

43.71 

43.7105 

-217.79 

-217.794 

62.642 

62.6422 

230.8 

230.8 

2 

68.795 

68.7948 

-104.51 

-104.509 

33.481 

33.4808 

129.52 

129.522 

3 

179.7 

179.699 

-81.806 

-81.8057 

72.166 

72.166 

210.22 

210.218 

4 

105.08 

105.078 

-140.79 

-140.792 

137.32 

137.316 

222.98 

222.979 

5 

11.764 

11.7636 

-125.14 

-125.139 

64.311 

64.311 

141.19 

141.188 

6 

-7.2061 

-7.20607 

-21.365 

-21.3654 

2.6549 

2.65487 

22.704 

22.7036 

7 

93.71 

93.7105 

-43.192 

-43.1922 

70.025 

70.0253 

124.7 

124.703 

8 

6.5107 

6.51075 

-35.082 

-35.0822 

158.41 

158.407 

162.38 

162.376 

9 

62.168 

62.1681 

-95.501 

-95.5015 

76.254 

76.2537 

137.11 

137.113 

Table  5.  Flux  Solution  Benchmark  Comparison  for  Thermal  Analysis 

(tetrahedral  element) 


7.2.2.2  Static  Analysis  Results: 

The  boundary  conditions  applied  to  the  model  for  the  static  analysis  case  consist  of  the  ba¬ 
sic  rigid  boundary  constraints.  All  three  degrees  of  freedom  for  nodel  are  constrained.  Node  2  is 
constrained  in  Y  and  Z  while  nodes  3  and  4  are  constrained  in  Z  only.  Additionally  static  load  of 
100  is  applied  to  nodes  5,  6,  7,  8  in  the  negative  Z  direction.  Table  6  and  Table  7  show  the  com¬ 
parison  of  the  displacement  and  stress  solutions  respectively  for  the  static  analysis. 
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NODE 

Ux 

U. 

ANSYS 

FEECAP 

ANSYS 

FEECAP 

ANSYS 

FEECAP 

1 

0 

0 

0 

0 

0 

0 

2 

8.667E-07 

8.667E-07 

0 

0 

0 

0 

3 

7.408E-07 

7.408E-07 

8.667E-07 

8.667E-07 

0 

0 

4 

-1.26E-07 

-1.26E-07 

8.667E-07 

8.667E-07 

0 

0 

5 

-8.31E-09 

-8.31E-09 

-8.31E-09 

-8.31E-09 

-4.33E-06 

-4.33E-06 

6 

6.945E-07 

6.945E-07 

1 .722E-07 

1.722E-07 

-3.39E-06 

-3.39E-06 

7 

7.491E-07 

7.491E-07 

8.75E-07 

8.75E-07 

-4.33E-06 

-4.33E-06 

8 

4.633E-08 

4.633E-08 

6.945E-07 

6.945E-07 

-3.39E-06 

-3.39E-06 

9 

3.704E-07 

3.704E-07 

4.333E-07 

4.333E-07 

-I.92E-06 

-1.92E-06 

Table  6.  Displacement  Solution  Benchmark  Comparison  for  Static  Analysis 

(tetrahedral  element) 


NODE 

CJ, 

<^z 

'^xz 

ANSYS 

FEECAP 

ANSYS 

FEECAP 

ANSYS 

FEECAP 

ANSYS 

FEECAP 

ANSYS 

FEECAP 

ANSYS 

FEECAP 

1 

1,8949 

1.8949 

1.8949 

1.8949 

113.42 

113.42 

1.3316 

1.3316 

-1.7387 

-1.7387 

-1,7387 

-1,7387 

2 

0.24033 

0.24033 

0.24033 

0.24033 

104.2 

104.2 

0.66581 

0.66581 

-0.69525 

-0.69525 

0.69525 

0.69525 

3 

1.8949 

1.8949 

1.8949 

1.8949 

113.42 

113.42 

1.3316 

1.3316 

1.7387 

1.7387 

1.7387 

1.7387 

4 

0.24033 

0.24033 

0.24033 

0.24033 

104.2 

104.2 

0.66581 

0.66581 

0.69525 

0.69525 

-0.69525 

-0.69525 

5 

0.16484 

0.16484 

0,16484 

0.16484 

99.311 

99.311 

-0.27858 

-0.27858 

-6.2989 

-6.2989 

-6.2989 

-6.2989 

6 

-1.9877 

-1.9877 

-1.9877 

-1.9877 

84.487 

84.487 

-1.6492 

-1.6492 

0.3071 

0.3071 

-0.3071 

-0.3071 

7 

0.16484 

0.16484 

0.16484 

0.16484 

99.3 1 1 

99.311 

-0.27858 

-0.27858 

6.2989 

6.2989 

6.2989 

6.2989 

8 

-1.9877 

-1.9877 

-1.9877 

-1.9877 

84.487 

84.487 

-1.6492 

-1.6492 

-0.3071 

-0.3071 

0.3071 

0.3071 

9 

0 

1.14E-14 

0 

-2.3E-15 

100 

100 

0 

5.69E-15 

0 

-3.33E-15 

0 

8.448E-15 

Table  7.  Stress  Solution  Benchmark  Comparison  for  Static  Analysis 

(tetrahedral  element) 


7.2.2.3  Combined  Thermal  Stress  Analysis  Results: 

The  combined  analysis  involves  imposing  thermal  and  static  boundary  conditions  on  the 
model.  To  test  the  robustness  of  the  formulation,  the  model  was  subjected  to  flux  boundary 
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condition  on  two  faces,  convection  on  two  other  faces,  internal  heat  generation  in  three  elements 


and  prescribed  temperatures  on  three  nodes  and  even  thermal  load  data  as  thermal  boundary  con¬ 
ditions.  As  static  boundary  conditions,  some  degrees  of  freedom  were  constrained  to  zero  dis¬ 
placement,  two  nodes  had  one  degree  of  freedom  each  with  prescribed  displacements  greater  than 
zero.  Tables  8  through  1 1  show  the  comparison  of  the  thermal  and  static  stress  solutions  obtained 
by  FEECAP  with  respect  to  ANSYS.  As  can  be  seen  from  the  temperature,  flux,  displacement  and 
stress  solutions,  there  is  a  good  match  of  results  between  ANSYS  and  FEECAP  2.5. 


NODE 

TEMPERATURE 

ANSYS 

FEECAP 

1 

78.514 

78.596 

2 

89.465 

89.332 

3 

76.153 

75.313 

4 

40 

40 

5 

86.076 

86.455 

6 

80 

80 

7 

91.295 

91.605 

8 

90 

90 

9 

77.415 

77.524 

Table  8.  Temperature  Solution  Benchmark  Comparison  for  Combined  Thermal  Stress  Analysis 

(tetrahedral  element) 
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NODE 

F 

LUX-X  F 

LUX-Y 

FLUX-Z 

TOTAL  FLUX 

ANSYS 

FEECAP 

ANSYS 

FEECAP 

ANSYS 

FEECAP 

ANSYS 

FEECAP 

1 

84.427 

83.445 

-162.9 

-163.216 

67.486 

69.057 

195.49 

195.886 

2 

86.53 

85.091 

-79.385 

-80.422 

45.424 

47.4331 

125.91 

126.326 

3 

163.72 

160.857 

-83.599 

-85.8041 

92.753 

97.1698 

205.91 

206.59 

4 

118.38 

116.796 

-111.24 

-112.126 

144.53 

146.32 

217.44 

218.227 

5 

56.307 

55.0152 

-72.086 

-72.9816 

70.83 

71.3372 

115.69 

115.94 

6 

18.84 

18.2317 

7.2556 

7.52225 

14.261 

14.0482 

24.717 

24.2143 

7 

90.601 

90.3017 

-37.792 

-37.695 

85.99 

88.205 

130.5 

131.74 

8 

31.586 

31.6876 

-5.4912 

-5.9336 

162.92 

162.379 

166.05 

165.548 

9 

80.839 

79.7359 

-67.377 

-68.0383 

86.48 

87.9061 

136.21 

136.801 

Table  9:  Flux  Solution  Benchmark  Comparison  for  Combined  Thermal  Stress  Analysis 

(tetrahedral  element) 


NODE 

1 

4 

Uy 

ANSYS 

FEECAP 

ANSYS 

FEECAP 

ANSYS 

FEECAP 

1 

0 

0 

0 

0 

0 

0 

2 

-1.08E-03 

-1.08E-03 

O.OOE+00 

O.OOE+00 

0 

0 

3 

-1.21E-03 

-1.21E-03 

-7.96E-04 

-7.94E-04 

0 

0 

4 

-2.02E-04 

-2.03E-04 

-1.05E-03 

-1.05E-03 

0 

0 

5 

3.88E-04 

3.88E-04 

6.66E-04 

6.66E-04 

lE-05 

lE-05 

6 

-5.88E-04 

-5.89E-04 

8.93E-04 

8.93E-04 

-8.69E-04 

-8.69E-04 

7 

-9.38E-04 

-9.40E-04 

2.00E-05 

2.00E-05 

-1.28E-03 

-1.28E-03 

8 

2.45E-04 

2.43E-04 

-3.72E-04 

-3.73E-04 

-9.00E-04 

-9.01E-04 

9 

-4.51E-04 

-4.53E-04  I 

-1.97E-04 

-1.97E-04 

-3.58E-04 

-3.58E-04 

Table  10.  Displacement  Solution  Benchmark  Comparison  for  Combined  Thermal  Stress  Analysis 

(tetrahedral  element) 
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NODE 

a, 

C 

'xy 

^yz 

<^xz 

ANSYS 

FEECAP 

ANSYS 

FEECAP 

ANSYS 

FEECAP 

ANSYS 

FEECAP 

ANSYS 

FEECAP 

ANSYS 

FEECAP 

1 

-4578.1 

-22722 

-22764 

2316.4 

-4652.7 

-4657 

-1590.2 

-1590.9 

2 

-10507 

-14999 

-132.53 

-5010.3 

-5016 

-758.73 

-751.93 

3 

-5053,7 

2181 

-1565.3 

-6366.6 

-6356.6 

-1368.2 

-1353.6 

4 

17951 

15349 

508,05 

-3900.7 

-3897.6 

-900.39 

-901.91 

5 

-7040 

-17669 

1066.2 

-1202.7 

-1204.5 

2067.8 

2066.5 

6 

-2116 

-2166,5 

-564.83 

-3207.1 

-3220.2 

1634.5 

1632.9 

7 

-13918 

-6805.5 

-1241.1 

-4814.1 

-4816.3 

-382.15 

-383.19 

8 

-9753.8 

-8023.3 

-343.29 

-1086.5 

-1086.7 

965,89 

952.07 

9 

-3497.7 

-5830.2 

-5884.8 

0 

-0.020847 

-3716.2 

-3718.2 

0 

l.OlE-04 

Table  1 1 .  Stress  Solution  Benchmark  Comparison  for  Combined  Thermal  Stress  Analysis 


(tetrahedral  element) 

In  summary,  all  results  for  temperatures,  fluxes,  displacements  and  stresses  agree  with  those 
obtained  from  FEECAP.  The  results  from  the  thermal  and  static  cases  agree  exactly  with  the  AN- 
SYS  results.  The  results  for  the  combined  thermal  and  static  analysis  also  closely  match  the  re¬ 
sults  from  ANSYS  with  a  maximum  error  of  1%  in  the  temperature  and  displacement  solution,  a 
maximum  of  4.8%  error  in  the  flux  solution  and  a  maximum  of  4%  error  in  the  stress  solution. 


8  Conclusions 

To  compare  the  relative  merits  of  various  alternative  MCM  package  designs,  it  is  necessary 
to  accurately  predict  the  thermal  behavior  of  small  integrated  circuit  die  features,  such  as  field  ef¬ 
fect  transitions,  and  to  accurately  predict  the  mechanical  behavior  of  small,  first-level  chip-to- 
substrate  interconnects.  Through  a  combination  of  finite  element  analyses,  submodeling  tech¬ 
niques,  and  the  application  of  specific  artificial  intelligence  techniques,  methodologies  were  de¬ 
veloped  and  implemented  in  a  proof-of-concept  software  tool,  called  the  Intelligent  Multichip 
Module  Analyzer,  to  predicting  Methodologies  were  developed  in  this  The  objective  of  this 


53 


project  was  to  develop  rapid,  accurate  methodologies  for  prediction  of  the  thermal  behavior  of  in¬ 
tegrated  circuit  features  and  mechanical  behavior  of  interconnects  and  to  implement  these  meth¬ 
odologies  into  the  Intelligent  Multichip  Module  Analyzer  (IMCMA).  The  combination  of  a 
high-level,  object-oriented  data  representation,  the  blackboard-based  problem  solving  paradigm, 
and  these  submodeling  methodologies  in  IMCMA  creates  a  powerful  prototype  computer-based 
tool  for  comparison  of  MCM  designs. 

IMCMA  uses  the  thermal  submodeling  process  to  rapidly  provide  an  MCM  designer  with 
an  estimate  of  the  maximum  package  temperature  of  a  device.  By  "zooming"  onto  critical  regions 
from  a  macroscopic  analysis  of  an  MCM  package,  IMCMA  avoids  the  time-consuming  process  of 
analyzing  large  finite  element  models.  With  IMCMA,  the  designer  can  investigate  the  thermal 
characteristics  of  several  candidate  MCM  designs  in  the  time  it  would  normally  take  to  analyze  a 
single  package. 

The  intercormect  submodeling  process  incorporated  in  IMCMA  also  uses  an  initial  mini¬ 
mum  degree  of  freedom  analysis  of  an  MCM  model  to  identify  a  critical  region.  Through  integra¬ 
tion  with  Pro/ENGINEER  and  the  addition  of  the  tetrahedral  element  to  the  finite  element  solver 
FEECAP,  IMCMA  has  the  capability  to  model  and  analyze  complex  interconnect  geometry.  Af¬ 
ter  the  creation  of  a  parametric  interconnect  model,  the  designer  can  complete  several  IMCMA 
runs  to  quickly  identify  an  MCM  configuration  that  produces  acceptable  stress  levels  in  its  inter¬ 
connects.  The  complex  geometry  of  MCM  interconnects  will  necessitate  large  finite  element 
models,  but  through  automatic  determination  of  the  critical  interconnect  and  its  boundary  condi¬ 
tions,  IMCMA  will  still  provide  substantial  time  savings  in  the  design  process. 
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Thermal  and  interconnect  submodeling  have  both  been  successfully  implemented  into  IM- 
CMA  to  decrease  the  amount  of  analysis  time  needed  to  complete  the  MCM  design  process.  The 
fact  that  these  submodeling  methodologies  are  automated  provides  time  savings  in  addition  to 
those  previously  discussed.  The  IMCMA  user  only  needs  to  provide  input  files  that  contain  high 
level  geometric  and  boundary  condition  information  to  run  IMCMA.  The  MCM  design  engineer 
doesn't  have  to  perform  the  tedious,  time-consuming  tasks  involved  in  the  creation  of  finite  ele¬ 
ment  models.  The  engineer's  time  can  be  used  more  productively  in  the  consideration  of  other  de¬ 
sign  concerns. 

9  Bibliography 

{\)  Proceedings  of  the  1992  IEEE  Multi-chip  Module  Conference,  March  18-20,  1992,  Santa 
Cruz,  CA,  p.  V. 

(2)  R.  R.  Tummala,  E.  J.  Rymaszewski,  editors.  Microelectronics  Packaging  Handbook,  Van 
Nostrand  Reinhold,  New  York,  NY,  1989,  pp.  391-409. 

(3)  K.  Azar,  "Thermal  Design  Considerations  for  Multichip  Module  Applications,"  Multichip 
Module  Technologies  and  Alternatives:  The  Basics,  D.  A.  Doane,  P.  D.  Franzon,  editors.  Van 
Nostrand  Reinhold,  New  York,  NY,  1993,  pp.  569-613. 

(4)  R.  R.  Brown,  J.  C.  Hodge,  "Engineering  Data  Management  Issues  and  Alternative 
Strategies,"  Managing  Engineering  Data:  The  Competitive  Edge,  edited  by  R.  E.  Fulton, 
1987  ASME  International  Computers  in  Engineering  Conference  and  Exhibition,  August 
9-13,  1987,  NY,  NY,  pp.  23-26. 

(5)  K.  H.  Law,  "Conceptual  Database  Design  for  Engineering  Modeling,"  Managing  Engineering 
Data:  Emerging  Issues,  1988  ASME  International  Computers  in  Engineering  Conference  and 
Exhibition,  July  31  -  August  4,  1987,  San  Francisco,  CA,  pp.  19-26. 


55 


(6)  R.  G.  G.  Cattell,  Object  Data  Management:  object-oriented  and  extended  database  systems, 
Addison-Wesley,  1992. 

(7)  D.  Corkill,  "Blackboard  ’Systems,"  AI  Expert,  Sept.  1991,  pp.  41-47. 

(8)  D.  Corkill,  "An  Architecture  for  Intelligent  Multichip  Module  Reliability  Analysis,"  a 
proposal  for  Rome  Laboratory,  NY. 

(9)  D.P.  Finn,  J.  B.  Crimson,  N.  M.  Harty,  "An  Intelligent  Mathematical  Modelling  Assistant  for 
Analysis  of  Physical  Systems,"  Proceedings  of  the  1992  ASME  International  Computers  in 
Engineering  Conference  and  Exhibition,  August  2-6,  San  Francisco,  CA,  Volume  2,  pp. 
69-74. 

(10)  Pepper,  D.W  and  Heinrich,  J.C,  The  Finite  Element  Method:  Basic  Concepts  and 

Hemisphere  Publishing  Corporation,  1992,  pp.  153-177. 

(11)  Zienkiewicz,O.C  and  Taylor,  R.L,  The  Finite  Element  Method,  Fourth  Edition,  vol.  1,  pp. 
89-95. 

(12)  Cook,R.D,  Malkus,D.S  and  Plesha,  M.E,  Concepts  and  Applications  Of  Finite  Element 
Analysis,  Third  Edition,  pp.  147-155. 

(13)  Shephard,  M.S.,  "Finite  Element  Grid  Optimization-  A  Review",  Finite  Element  Grid 
Optimization,  M.S. Shephard  and  R.Gallagher.  eds.,  ASME  Special  Publication  PVP-38, 
N.Y.,  1979. 

(14)  Babuska,  I.,  O.C.  Zienkiewicz,  J.  Gago,  and  E.R.A.  Oliveira,  Accuracy  Estimates  and 
Adaptive  Refinements  in  Finite  Element  Computations,  John  Wiley  &  Sons,  New  York,  1986. 

(15)  Ahmed  K.  Noor,  and  Ivo  Babuska,  "Quality  Assessment  and  Control  of  Finite  Element 
Solutions,"  International  Journal  of  Applied  Finite  Elements  and  Computer  Aided 
Engineering,  Vol.  3,  No.  1,  April  1987. 


56 


(16)  Babuska,  I.,  and  W.C.  Rheinboldt,  "A  Posteriori  Error  Estimations  for  The  Finite  Element 
Method,"  International  Journal  for  Numerical  Methods  in  Engineering,  Vol.  12, 
pp.1597-1615,  1978. 

(17)  Rheinboldt,  W.C.,  "Error  Estimates  for  Nonlinear  Finite  Element  Computations", 
Computers  &  Structures,  Vol.  20,  No.  1-3,  pp.  91-98,  1985. 

(18)  Babuska,  I.  and  A.  Miller,  "The  Post-processing  Approach  in  the  Finite  Element  Method  - 
Part  3;  A  Posteriori  Error  Estimates  and  Adaptive  Mesh  Selection",  International  Journal  for 
Numerical  Methods  in  Engineering,  Vol.  20,  pp.231 1—2324,  1984. 

(19)  Kelly,  D.W.,  J.P.  Gago,  O.C.  Zienkiewicz,  and  I.  Babuska,  "A  Posteriori  Error  Analysis  and 
Adaptive  Processes  in  the  Finite  Element  Method;  Part  I,  Error  Analysis",  International 
Journal  for  Numerical  Methods  in  Engineering,  Vol.  19,  pp.  1593-1619,  1983. 

(20)  Gago,  J.P.,  D.W.  Kelly,  O.C.  Zienkiewicz,  and  I.  Babuska,  "A  Posteriori  Error  Analysis  and 
Adaptive  Processes  in  the  Finite  Element  Method:  Part  II,  Adaptive  Mesh  Refinement", 
International  Journal  for  Numerical  Methods  in  Engineering,  Wo\.  19,  pp.  1621—1656,  1983. 

(21)  Ainsworth,  M.,  J.Z.  Xhu,  A.W.  Craig  and  O.C.  Zienkiewicz,  "Analysis  of  the 
Zienkiewicz-Zhu  a  posteriori  error  estimator  in  the  finite  element  method".  International 
Journal  for  Numerical  Methods  in  Engineering,  28,  pp.  2161-2174,  1989. 

(22)  Zienkiewicz,  O.C.  and  Z.  Zhu,  "A  Simple  Error  Estimator  and  Adaptive  Procedure  for 
Practical  Engineering  Analysis",  International  Journal  for  Numerical  Methods  in 
Engineering,  Vol.  24,  pp. 337— 357,  1987. 

(23)  Huang,  H.C.  and  R.W.  Lewis,  "Adaptive  Analysis  for  Heat  Flow  Problems  using  Error 
Estimation  Techniques",  Sixth  International  Conference  for  Numerical  Methods  in  Thermal 
Problems,  U.K.,  July  1989,  Pineridge  Press,  Swansea,  U.K.,  1989. 

(24)  Lewis,  R.W.,  H.C.  Huang,  A.S.  Usmani  and  J.T.  Cross,  "Finite  Element  Analysis  of  Heat 
Transfer  and  Flow  Problems  using  Adaptive  Remeshing  including  Application  to 


57 


Solidification  Problems",  International  journal  for  Numerical  Methods  in  Engineering,  Vol. 
32,  pp.767-781,  1991. 

(25)  Grosse,  I.R.,  P.  Katragadda  and  J.  Benoit,  "An  Adaptive  Accuracy  Based  A  Posteriori  Error 
Estimator",  Finite  Elements  in  Analysis  and  Design,  pp  70-95,  September,  1992. 

(26)  Bathe  and  Wilson,  "Numerical  Methods  in  Finite  Element  Analysis,"  Prentice-Hall,  New 
Jersey,  1976,  Chap.  6 

(27)  Zienkiewicz,  O.C.,  The  Finite  Method,  Mcgraw-Hill,  New  York,  3rd  edition,  1977 

(28)  Zienkiewicz,  O.C.  and  Z.  Zhu,  "Adaptivity  and  Mesh  Generation",  International  Journal  for 
Numerical  Methods  in  engineering,  Vol.  32,  pp.  783-810,  1991. 

(29)  Blacker,  T.D.  and  M.B.  Stephenson,  "Paving:  A  New  Approach  to  Automated  Quadrilateral 
Mesh  Generation,"  Sandia  Report,  SAND90-0249.  UC-705,  1990. 

(30)  Blacker,  T.D.,  J.  Jung  and  W.R.  Witkowski,  "An  Adaptive  Finite  Element  Technique  Using 
Element  Equilibrium  and  Paving,"  ASME  Paper  NO.  90-WA/CIE-2 

(31)  Grosse,  I.R,  Katragadda,  P  and  Jog,  A,  "Knowledge  Sources  For  An  Intelligent  MCM 
Analyzer,"  Final  Technical  Report,  Rome  Laboratory  Airforce  Materiel  Command,  Griffiss 
Air  Force  Base,  New  York,  June  1993 

(32) Doane,  D.  A.  and  Franzon,  P.  D.,  editors.  Multichip  Module  Technologies  and  Alternatives: 
The  Basics,  Van  Nostrand  Reinhold,  New  York,  NY,  1993 


58 


Appendix  A:  A  Posteriori  Error  Estimation  and  Adaptive  Mesh 
Refinement  for  Combined  Thermal-Stress  Finite  Element  Analysis 

Previous  sections  have  dealt  with  how  the  critical  region  of  a  model  can  be  identified  by 
performing  a  macro  analysis  of  the  complete  model  and  then  doing  an  automatic  modeling  and 
analysis  of  submodels.  It  is  to  be  noted  that  there  exists  inherent  discretization  errors  in  the  FEA 
method  and  hence  it  is  necessary  for  the  user  to  be  aware  of  the  magnitude  of  such  errors  and  their 
effect  on  the  solution.  Also,  in  the  domain  of  MCMs,  all  the  components  are  thermally  loaded  and 
statically  constrained  and  hence  there  is  a  need  for  an  error  estimator  which  can  compute  the  mag¬ 
nitude  of  errors  for  each  element  under  such  type  of  dual-loading. 

A,  1.  Literature  review 

In  the  past  few  years  a  significant  amount  of  research  has  been  done  in  the  area  of  a  priori 
and  a  posteriori  error  estimations  for  the  finite  element  solution.  Various  methods  have  been  pro¬ 
posed  and  some  of  them  have  been  tested  rigorously.  For  a  more  complete  review  of  related  work 
done  so  far  in  error  estimations  and  adaptive  remeshing  the  reader  is  referred  to  [13]-[15].  A  ma¬ 
jority  of  the  algorithms  developed  have  been  for  static  stress  analysis. 

Much  of  the  earlier  work  in  residual  based  error  analysis  was  predominantly  mathematical 
and  was  done  by  Babuska  and  associates  [16]-[18].  Based  on  the  residual  of  the  differential  equa¬ 
tion,  the  authors  derived  error  bounds  for  the  energy  norm  of  the  error  and  element  error  indica¬ 
tors  were  introduced  to  determine  which  elements  were  to  be  refined.  Based  on  this  work  Kelly  et 
al  [19]  proposed  a  residual  based  error  estimator  which  used  special  hierarchical  shape  functions. 
These  hierarchical  shape  functions  permitted  efficient  p-based  error  estimates.  A  subsequent  paper 
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by  Gago  et  al  [20]  was  published  which  offered  several  strategies  for  using  an  error  estimate  to  re¬ 
fine  a  mesh. 

Subsequently,  Zienkiewicz  and  Zhu  outlined  a  simple  stress-based  error  estimator  and  asso¬ 
ciated  adaptivity  algorithm.  This  method  involves  obtaining  a  global  least-squares  fit  of  the  dis¬ 
continuous  finite  element  stress  field  with  a  C“  continuous  stress  field.  This  C°  continuous  stress 
field  is  taken  as  an  approximation  to  the  true  stress  field  and  the  difference  between  the  two  stress 
fields,  as  measured  by  the  energy  or  L2  norm,  represents  an  estimate  of  the  discretization  error. 

Huang  and  Lewis  [23]  developed  an  error  estimation  and  adaptive  technique  for  linear 
steady-state  conduction  problems.  In  a  subsequent  paper  Lewis  et  al  [24]  extended  the  algorithm 
to  the  solution  of  non-linear  transient  heat  conduction  problems.  For  a  measure  of  the  discretiza¬ 
tion  errors,  the  authors  used  heat  flux  which  is  a  function  of  temperature  gradients.  For  calculating 
the  error  norms,  they  adopted  the  idea  presented  by  Zienkiewicz  and  Zhu  [22].  Zienkiewicz  and 
Zhu  [22]  showed  that  globally  smoothed  values  of  stresses  are  nothing  but  higher  order  approxi¬ 
mation  to  the  stresses  obtained  from  a  finite  element  analysis,  can  be  used  instead  of  the  exact 
stresses.  The  requirement  that  for  a  near  optimal  mesh  all  the  elements  of  the  final  mesh  must 
contain  an  approximately  equal  error  was  used. 

Grosse  et  al  [25]  have  proposed  a  stress  based  finite  element  error  estimator  with  a  new  con¬ 
cept  of  adaptive  accuracy  and  mesh  optimality.  They  have  shown  that  the  meshes  obtained  from 
the  /7-based  adaptivity  algorithm  were  more  accurate  in  the  highly  stressed  regions  of  the  domain 
with  substantially  reduced  degrees  of  freedom  compared  to  the  adaptive  meshes  obtained  with  an 
uniform  accuracy  criteria.  This  was  achieved  by  formulating  the  accuracy  requirement  to  be  a 
function  of  the  relative  behavior  of  the  stress  field  distribution  across  the  domain.  That  is,  the  net 
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result  of  this  adaptive  accuracy  criteria  is  a  nonuniform  distribution  of  discretization  error  with 
minimum  error  in  maximum  stress  regions  and  visa  versa. 

Here  we  are  concerned  with  coupled  thermal  stress  analysis  problems.  For  such  dual¬ 
loading  type  problems  a  thermal  analysis  is  first  conducted,  followed  by  a  static  stress  analysis. 
The  two  analysis  are  coupled  to  each  other  via  the  thermal  strains  from  the  thermal  analysis  which 
play  the  role  of  initial  strains  in  the  static  analysis.  Thus,  the  two  physical  phenomenon  are  unidi- 
rectionally  eoupled.  That  is,  although  thermal  behavior  directly  affects  elastostatic  behavior  via 
thermal  strains,  the  elastostatic  behavior  is  considered  to  have  a  negligible  effect  on  the  thermal 
behavior.  This  is  the  case  for  most  thermal-stress  engineering  problems,  but  bi-directional 
thermal-elastostatic  coupling  can  exist  for  special  problems,  especially  problems  involving  ther¬ 
mal  contact  resistance. 

When  performing  the  error  analysis,  the  error  algorithm  can  be  applied  at  two  different 
stages.  It  can  be  applied  to  the  thermal  analysis  results  or  it  can  be  conducted  after  the  statie  stress 
analysis  is  completed.  It  is  important  to  note  that  for  a  given  finite  element  mesh,  the  thermal 
analysis  solution  and  the  stress  analysis  solution  will  yield,  in  general,  different  diseretization  er¬ 
rors  throughout  the  domain.  Furthermore,  since  the  stress  solution  is  coupled  to  the  thermal  solu¬ 
tion  (via  thermal  strains)  the  thermal  solution  discretization  errors  can  adversely  affect  the 
elastostatic  solution.  This  effect  will  be  undetected  in  a  posteriori  error  analysis  of  the  elastostatic 
solution. 

The  most  rigorous  solution  to  this  problem  is  to  first  carry  out  an  iterative  a  posteriori  error 
analysis  and  adaptive  mesh  refinement  process  for  the  thermal  analysis  before  proceeding  to  the 
stress  analysis.  Then  an  iterative  a  posteriori  error  analysis  and  adaptive  mesh  refinement  is 
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conducted  for  the  stress  analysis  portion  of  the  problem.  This  strategy,  however,  is  prohibitively 
computationally  expensive.  Clearly  two  separate  analysis,  error  estimation  and  adaptive  mesh  re¬ 
finement  iteration  loops  are  required.  Moreover,  for  each  elastostatic  analysis  (with  the  exception 
of  the  first  one)  elemental  thermal  strains  computed  from  the  thermally  converged  mesh  solution 
must  be  mapped  onto  the  current  elastostatic  mesh. 

Simpler,  more  computationally  efficient  schemes  are  clearly  needed.  In  this  paper  we  pre¬ 
sent  error  estimators  for  combined  thermal  and  elastostatic  analysis  and  an  adaptive  remeshing 
methodology  which  effectively  refines  the  mesh  based  on  both  error  estimators  in  a  cost  effective 
manner. 


A.2.  Error  estimator  and  adaptive  accuracy  criteria  for  static  stress  analysis: 

In  previous  work  [25]  we  introduced  a  simple  error  estimator  for  static  analysis,  similar  to 
the  Zienkiewicz  -Zhu  error  estimator  [21],  based  on  the  effective  stress  (von  Mises)  function 
and  the  Lj  error  norm. 


I|£(‘<^v)||2  =  || 


[(a.- 


(17) 


where  is  the  effective  stress  function  directly  obtained  from  the  finite  element  displacement 
and  strain  solution  and  ‘’a*  is  a  C“  continuous  'improved'  effective  stress  function  computed  via 

stress  recovery  techniques  presented  in  reference  [25]. 

The  relationship  of  the  effective  stress  function  to  both  distortional  strain  energy  and  to  the 
distortional  energy  failure  theory  is  well  known.  Also,  it  has  been  shown  that  the  L,  norm  of  the 
element's  von  Mises  stress  field  is  proportional  to  the  square  root  of  the  element's  distortional 
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strain  energy  [25].  Hence,  a  prudent  choice  is  to  base  the  error  estimation  algorithm  on  the  ele¬ 
ment's  effective  stress  function  and  use  it  as  a  measure  of  the  element's  discretization  error. 

The  elemental  error  in  the  norm  sense  is  forced  to  be  less  than  some  small  fraction  of  a  ref¬ 
erence  L2  norm  value: 

<12  <  r|(||av|l2);;^y  (18) 

where  r\  is  small.  If  the  reference  value  is  a  local  (i.e.  element  based)  value,  then  excessive  refine¬ 
ment  due  to  a  constant  relative  accuracy  requirement  will  be  imposed  for  all  elements  regardless 
of  the  stress  state  of  the  element.  On  the  positive  side,  an  element  based  reference  value  for  the  Lj 
norm  bounds  the  global  L2  error  norm  in  to  be  small  compared  to  the  global  Lj  norm  in  a^.  On 
the  other  hand,  if  a  single  global  value  is  taken  as  the  reference,  then  excessive  mesh  refinement  is 
avoided,  but  the  global  Lj  error  norm  in  may  not  be  small  (i.e..  bounded)  compared  to  the 
global  L2  norm  of  a^.  This  possibility  can  be  minimized  by  adopting  a  bifunctional  adaptive  refer¬ 
ence  norm  as  follows: 

(llavllz)*,/"  max  (11  •o;||z,  G/;Ms(||  c:||z)  )  (19) 

where  is  the  global  root  mean  square  of  all  the  elemental  norms  in  the  effective 

stress: 

=  [i  \Kf2dv\1A' 

_Ve=l  ey  /  J  Le=l  - 
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Note  that  Eq.  (19)  involves  a  combination  of  the  elemental  and  global  effective  stress  norms 
in  such  a  way  as  to  preserve  the  relative  error  norms  in  regions  of  interest,  namely,  in  elements  of 
high  stress  states,  while  eliminating  the  problem  of  excessive  mesh  refinement  in  regions  highly 
understressed.  Figure  16  illustrates  the  bi-functional  adaptive  reference  norm. 


Figure  16:  Illustration  of  the  bifunctional  reference  norm 

Substitution  of  Eq.  (19)  into  (18)  results  in  a  bi-functional  accuracy  requirement  based  on 

the  Lj  norm  of  the  elements  effective  stress  state  compared  to  the  global  root  mean  square  value 

of  all  the  element  L,  effective  stress  norms: 

|hav-^a:||<Ti  max(pa:||2,Gi?Ms(||a:||2))  (21) 

Eq.  (21)  constitutes  a  bi-functional  adaptive  accuracy  requirement.  In  the  stress  based  error 

estimator  and  adaptivity  algorithms  presented  by  Grosse  et  al  [25]  introduced  a  new  concept  of 

adaptive  accuracy  based  on  the  idea  that  an  optimal  mesh  from  an  engineering  viewpoint  is  one 

which  adapts  according  to  the  element's  stress  state.  This  means  that  the  desired  accuracy  is 
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achieved  at  regions  of  interest  with  the  least  number  of  degrees  of  freedom.  Typically  this  corre¬ 
sponds  to  tighter  accuracy  requirements  in  regions  of  higher  stress  states  and  relaxed  accuracy  re¬ 
quirements  in  regions  of  lower  strqss  states. 

An  adaptive  accuracy  can  be  extended  to  the  element  level  by  allowing  r\,  the  target  relative 
error  norm  ratio  to  adapt  to  the  element  stress  state  as  follows,  by  introducing  an  elemental  nomi¬ 
nal  target  accuracy  based  on  the  element  stress 


GRMS(gI) 

RMSi 


(22) 


In  the  above  expression  r\  is  the  nominal  target  accuracy,  GRJVIS(g/  )  is  the  global  RMS  measure 
of  the  effective  stress  function  over  the  entire  domain  and  RMS(‘^ct^,*  )  is  the  RMS  of  the  effective 
stress  function  distribution  for  element  e: 


GRMSiGl) 


RMSi^Gl) 


l(G;fdV/V 


I  \(Giydv\iv 

.e=l  cy 


15 


j  yG,fdv\i‘v 

\  cy 


_  I 


(23) 

(24) 


Combining  the  bi-functional  adaptive  reference  norm  with  the  elemental  adaptive  target  ac¬ 
curacy,  the  accuracy  criteria  becomes 

max  (\\  '"(7*112,  GRMs{\\  aClQ  )  (25) 


< 


grms{g;) 
RMSGg:)  J 


Experimentally,  this  modification  has  proven  quite  promising,  limiting  the  unjustifiably 
high  values  for  error  measurement  parameters  encountered  with  the  element  reference  value  in 
understressed  regions. 
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In  Eq.  (22),  a  is  a  constant  >0  which  controls  the  degree  of  accuracy  adaptivity.  If  a  =  0, 
then  a  nonadaptive  accuracy  requirement  is  imposed.  It  is  observed  that  the  effect  of  Eq.  (22)  is  to 

increase  the  accuracy  requirement  for  elements  with  RMS(%^,*  )  greater  than  the  global  RMS 
measure  of  the  smoothed  effective  stress  function.  Conversely,  the  accuracy  requirement  is  de¬ 
creased  for  all  elements  where  RMS(M^,*  )  <  GRMS(<Ty’ ). 

Praetical  bounds  must  be  placed  on  '"ri  computed  by  Eq.  (22)  to  prevent  mesh  transitioning 
problems  in  an  h-refmement  meshing  scheme.  The  user  specifies  a  priori  and  ‘'rimax^  the  tight¬ 
est  and  slackest  accuraey  requirements,  which  correspond  to  the  maximum  and  minimum  elemen¬ 
tal  RMS(%y*  )  values.  Henee,  from  Eq.  (22)  the  accuracy  requirements  on  the  element  with  the 
highest  stressed  state  and  the  element  with  the  least  stressed  state  are 


and 


•rimin  =  ri 
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A  relationship  between  and  a  can  be  obtained  by  dividing  Eq.  (26)  by  Eq.  (27), 


^  min 
^Pmax 


(28) 


Solving  Eq.  (28)  for  a  unique  value  of  a 


a 


(29) 


which  is  based  on  the  user  specified  values  for  and  and  the  calculated  minimum  and 


maximum  elemental  RMS  values  in  effective  stress.  Once  a  is  calculated  via  Eq.  (29)  it  can  be 
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substituted  into  either  Eq.  (26)  or  (27)  to  determine  the  value  of  t].  From  Eq.  (18),  a  dimensionless 
parameter  called  the  element  error  ratio  for  the  static  stress  solution  can  be  defined  as; 


‘'Ov-  ''a*ll2 


Ti(llavll2)«,^ 

and  substituting  Eq.(22)  in  Eq.  (30),  the  final  expression  for  %is  given  as: 


‘'av- 


f  GRMS(a:.)]  “  II  II 

If  >  U  then  the  element  size  di  is  decreased  and  if  *=^5  <  E  the  element  size  "h  is  increased.  If 
=  1  for  all  elements,  then  the  adaptive  accuracy  criteria  has  been  met  everywhere. 

A.3.  Error  estimator  for  thermal  analysis: 

Analogous  to  the  choice  of  the  effective  stress  fianction  as  the  basis  for  measuring  the  dis¬ 
cretization  error  in  static  analysis,  the  thermal  heat  flux  magnitude  is  chosen  for  error  analysis  in 
the  thermal  solution.  Following  the  same  steps  as  shown  above,  the  error  ratio  for  thermal  analy¬ 


sis  is  defined  as 


^q-  ^q*\\2 


(  GRMSiq 


The  heat  flux  q  is  given  by: 


{q}=-[K]YT 


and  the  heat  flux  magnitude  is: 

q  =  [{(lV{(l}Y  (34) 

where  [K]  is  the  thermal  conductivity  matrix  of  the  material  and  Y  T  the  gradient  in  the  tempera¬ 
ture  field  is  given  by 


V  T  = 
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Just  as  in  static  analysis  where  mesh  refinement  is  concentrated  in  regions  having  high 
stress  gradients,  in  thermal  analysis  mesh  refinement  will  be  seen  to  be  concentrated  in  regions  of 
high  heat  flux  gradients. 

A.4.  Error  estimator  for  coupled  thermal  and  static  analysis: 

When  a  body  is  subjected  to  both  thermal  and  mechanical  loads,  the  stress-strain  relations 
must  include  the  effects  of  thermal  expansion  (and  contraction)  by  subtracting  the  thermal  strains 

{Sq}  from  the  total  strains  {s}.  Therefore, 

{a}  =[£]({s}-{£o})  (36) 

where 

{£o}=a.Ar[  1  1  1  0  0  0  (37) 

for  a  three  dimensional  problem.  Here  a  is  the  coefficient  of  thermal  expansion  for  an  isotropic 
material  and  AT  is  the  temperature  field  relative  to  a  reference  temperature  at  a  stress  free  state. 
The  total  strains  are  due  to  both  the  thermally  induced  and  stress  induced  strains.  In  the  standard 

finite  element  formulation  the  equivalent  static  loads  induced  by  thermal  strains  are: 

{^Rt}=  j  VBYVE\{z^}dV  (38) 

‘V 

where  is  the  element  strain  displacement  matrix  which  relates  the  element  strain  field  to 

nodal  displacements  and  is  given  by  the  gradient  of  the  shape  function  matrix  [TJ]: 

{^8}  =  (39) 

where  Y 

The  static  analysis  system  of  equations  becomes 

[Ks]{D}  ^  {Rs} +{Rt}  (40) 
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Solving  Eq.  (40)  and  using  Equations  (39)  and  (36)  we  obtain  the  stress  solution  which  includes 
the  effects  of  the  thermally  induced  initial  strains. 

A.  5.  Methods  of  Mesh  Refinement: 

The  most  rigorous  approach  (Method- 1)  is  to  carry  out  an  iterative  a  posteriori  error  analy¬ 
sis  and  adaptive  mesh  refinement  process  for  the  thermal  analysis  before  proceeding  to  the  stress 
analysis,  and  then  perform  the  iterative  a  posteriori  error  analysis  and  adaptive  mesh  refinement 
process  for  the  static  stress  analysis  portion  of  the  problem.  By  this  method  we  first  refine  the 
mesh  to  converge  to  a  thermal  solution  within  a  user  specified  accuracy  level  and  then  perform 
further  refinement  to  obtain  a  stress  solution  within  an  acceptable  level  of  accuracy.  Since  this 
strategy  is  computationally  expensive,  it  is  not  recommended. 

A  second  method  (Method-2)  is  to  apply  elastostatic  error  analysis  after  both  the  thermal 
and  stress  analysis  are  conducted.  This  method  is  not  recommended  due  to  the  coupling  of  the 
stress  solution  to  the  thermal  solution  (via  thermal  strains)  and  the  effect  of  the  thermal  discretiza¬ 
tion  errors  in  the  stress  solution  which  can  go  undetected  in  an  a  posteriori  error  analysis  of  the 
elastostatic  solution. 

Instead,  we  propose  here  a  simple  scheme  referred  to  as  the  Proposed  Method.  The  thermal 
analysis  is  first  conducted  and  an  a  posteriori  error  analysis  on  this  thermal  solution  is  performed. 
From  this  thermal  error  analysis  we  obtain  the  corresponding  thermal  analysis  error  ratios 
Next,  a  elastostatic  analysis  is  done  using  the  thermal  strains  as  initial  strains  and  an  independent 
a  posteriori  error  analysis  is  performed  on  the  static  solution  to  obtain  the  corresponding  static 
analysis  error  ratios  We  thus  obtain  two  sets  of  elemental  error  ratios,  one  from  the  error 
analysis  conducted  on  the  thermal  solution  and  the  other  from  the  error  analysis  conducted  on  the 
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static  stress  solution  These  two  error  ratios  for  each  element  signify  the  aeeuracy  achieved  in 
the  thermal  solution  and  the  accuracy  achieved  in  the  stress  solution. 

A  logical  and  conservative  approach  is  to  use  the  greater  of  the  two  error  ratios  for  a  given 
element  as  the  basis  for  mesh  refinement.  Thus,  if  an  element  is  less  accurate  in  its  thermal  solu¬ 
tion  as  compared  to  its  stress  solution,  the  element's  thermal  analysis  error  ratio  is  used  for  recom¬ 
puting  the  new  element  size.  Hence  the  final  elemental  error  ratio  (%c)  due  to  the  combined  effect 

of  thermal  and  mechanical  (static)  loads  is  given  by  the  relation: 

^^c  =  max(^^7-/^5)  (41) 

An  adaptive  mesh  refinement  is  then  performed  based  on  these  error  ratios  (^^c)  which  take 
into  account  the  errors  in  both  the  thermal  and  static  stress  analysis  solutions.  Intuitively,  we  ex¬ 
pect  a  faster  convergence  with  this  strategy. 

A.  6.  New  element  size: 


From  the  asymptotic  convergence  rate  criteria,  one  can  derive  the  following  relationship  be¬ 


tween  new  and  old  element  sizes  [29]: 


(42) 


where  "jh  denotes  the  element  size  for  the  i“'  mesh,  p  is  the  order  of  the  polynomials  used  in  the 


shape  function  expansions  for  the  temperature  and  displacement  fields  N  and  A,  is  the  intensity  of 
singularities  (if  present). 

A.  7.  Implementation: 

The  error  analysis  algorithms  were  implemented  into  a  FORTRAN-77  code  called  FEE- 
CAP  which  is  a  three-dimensional  Finite  Element  Error  Control  and  Analysis  Package  for  steady 
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state  thermal  analysis,  static  analysis  and  combined  thermal-stress/static  analysis.  The  error  analy¬ 
sis  algorithm  in  FEECAP  estimates  the  finite  element  discretization  errors  and  based  on  these  es¬ 
timates  and  user  specified  accuracy  requirements,  computes  the  element  parameters  called  error 
ratios.  These  error  ratios  in  conjunction  with  an  automatic  mesh  generator  can  be  used  for  auto¬ 
matic  adaptive  mesh  refinement.  Since  the  adaptivity  algorithms  have  not  yet  been  developed 
within  the  IMCMA  environment,  the  implementation  of  the  above  adaptive  remeshing  algorithms 
for  the  map  meshing  have  been  limited  to  a  2-D  automatic  mesh  generation  code  called  FASTQ 
developed  at  Sandia  National  Laboratories.  A  schematic  flowchart  of  FEECAP  is  shown  in  Fig¬ 
ure  17. 

To  study  the  convergence  of  the  global  norm  of  the  errors  in  the  thermal  flux  field  and 
the  effective  stress  field,  we  define  global  error  norms  as: 


The  effectivity  index  (denoted  by  El  or  0)  is  a  dimensionless  parameter  which  measures  the 


effectiveness  of  the  error  estimator  and  is  defined  as 


/7  ^  P  * 


‘'Gv-  ‘'g, 


where  %  is  the  "exact"  solution  for  the  element.  The  effectivity  index  (El)  can  either  be  used  to 


monitor  error  measurement  for  an  individual  element  of  the  model,  or  a  measure  of  the  effective¬ 


ness  of  the  error  estimator  over  the  entire  domain  is  given  by  the  global  effectivity  index: 


numc  ,,  ^  ^  ,,9 

s  rav- 

e=\ 


Z  “^av-  ‘"cTv 


2 


(46) 
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USER 


S  =  Static  Analysis  only 
C  ==  Combined  (thermal  and  static  )  analysis 


Fig.  17:  Flowchart  of  FEECAP 
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Eqs.  (45)  and  (46)  can  be  written  out  in  a  similar  manner  for  the  thermal  solution  simply  by  re¬ 
placing  a,,  with  the  thermal  flux  vector  {q}  or  the  scalar  flux  magnitude  function  q„,. 

A. 8.  An  Example  Problem: 

A  simple  example  of  a  chip  with  prescribed  thermal  loads  and  static  constraints  was  se¬ 
lected  to  demonstrate  the  convergence  and  effectiveness  of  the  error  algorithm  outlined  earlier. 
Figure  18  shows  the  example  chosen  with  the  boundary  conditions  applied.  In  IMCMA,  nine 
analysis  runs  were  performed  with  default-number-of-elements  set  to  1,2,3, •■■•and  9  for  each  of  the 
runs.  The  ninth  mesh  with  729  elements,  1512  active  thermal  degrees  of  freedom  and  4692  active 
static  degrees  of  freedom  was  found  to  be  sufficiently  accurate  to  "truly"  represent  the  thermal 
flux  and  stress  distributions.  Figures  (19)  and  (20)  show  the  convergence  of  the  global  error 
norms  (calculated  as  shown  in  Eqs.  (43)  and  (44))  versus  the  active  degrees  of  freedom. 

Figs.  (21)  and  (22)  show  the  convergence  of  the  global  effectivity  index,  computed  by  Eq. 
(46)  for  the  thermal  and  static  stress  solutions.  These  graphs  illustrate  the  decay  of  the  maximum 
and  minimum  effectivity  indices  across  the  whole  domain  of  the  chip  example  for  each  mesh. 
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Uniform  thickness  =  1.0 
Thermal  loads:  applied  flux,  q  =  0.03 
prescribed  T  =  30 
all  other  surfaces  are  insulated 

Static  constraints:  x,  y  and  z  constraints  on  4  comers  of  bottom  surface 


Figure  18:  Example  of  a  chip  with  applied  thermal  and  static  boundary  conditions 
used  for  benchmarking  the  combined  error  estimation  algorithm. 
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Figure  21 :  Convergence  of  effectivity  index  (El)  versus  degrees  of  freedom  (DOF) 
for  thermal  solution 
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Figure  22:  Convergence  of  effecti' 
for  the  stress  solution 
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index  (El)  versus  degrees  of  freedom  (DOF) 


Appendix  B.  Running  Thermal  Submodeling  in  IMCMA 

RTJN-/MCM4 

Five  keywords  have  been  implemented  in  the  run-imcma  function  to  carry  out  submodeling 
within  the  Intelligent  Multichip  Module  Analyzer  {IMCMA)  software  system.  These  keywords 
are  optional  arguments  to  the  run-imcma  function. 


thermal-submodel  This  keyword  is  used  to  run  thermal  submodeling  in  IMCMA.  This 

includes  the  macroscopic  model  analysis  and  two  submodel  analyses. 

Example:  (run-/mcma  "example"  :submodel  t) 

The  submodeling  process  will  be  completed  from  input  information 
found  in  example.lisp.  The  "t"  parameter  added  to  the  :submodel 
option  stands  for  "true". 

interconnect  This  keyword  is  used  to  identify  the  critical  interconnect  on  an  MCM, 

and  to  analyze  the  Pro/Engineer  mesh  of  the  interconnect.  The 
displacements  at  the  critical  interconnect  locations  are  applied  as 
boundary  conditions  to  the  mesh  of  the  interconnect. 


Example:  {x\xn-imcma  "example2"  :intercormect) 

The  critical  interconnect  will  be  identified  from  the  model  and 
interconnect  specifications  in  example2.1isp.  The  Pro/Engineer  mesh 
will  then  be  analyzed  and  the  results  imported  into  the  IMCMA 
database. 


die-attach-material  and  die-attach-thickness 

These  keywords  must  be  used  together  to  define  the  two  die  attach 
parameters.  The  die  attach  parameters  are  used  in  two  different  ways. 
If  the  equivalent-k  keyword  option  is  used  in  run-imcma,  the 
parameters  are  used  to  calculate  equivalent  thermal  conductivities.  If 
the  submodel  option  is  used,  the  die  attach  is  actually  modeled  in  the 
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first  submodel.  The  parameters  are  then  used  to  define  a  die  attach 
component  and  to  bind  a  material  to  this  component 

Example:  (run-/wcma  "examples"  :submodel  t 

:die-attach-material  ':epoxy 
;die-attach-thickness  0.05) 

An  epoxy  die  attach  with  a  thickness  of  0.05  will  be  included  in  the 
first  submodel  analysis.  The  material  must  be  defined  in  the  database 
or  an  error  will  result. 

equivalent-k  This  keyword  instructs  the  system  to  create  materials  with  equivalent 

thermal  conductivities  for  the  macroscopic  model  analysis.  The 
effects  of  the  chip  and  the  die  attach  are  represented  by  giving  the 
chip  a  material  with  an  equivalent  conductivity.  The 
die-attach-material  and  die-attach-thickness  keywords  must  be 
included  in  the  run-imcma  command  for  this  option  to  run. 

Example:  {mn-imcma  "example4"  :submodel  t  :equivalent-k  t 

:die-attach-thickness  0.05 
:die-attach-material  ':epoxy) 

Chip  materials  will  be  replaced  with  materials  having  a 
conductivityequivalent  to  the  chip/die-attach  combination  for  the 
macroscopic  analysis.  The  original  materials  are  returned  for  the 
subsequent  submodel  analyses,  and  a  die  attach  is  modeled  for  the 
first  submodel  analysis. 

COMPONENTS  AND  THERM  AT.  STIRMODRTJNG 

When  running  a  thermal  submodel  analysis,  both  components  and  2nd-subcomponents  must 
be  defined  by  using  the  defcomponent  and  defindsubcomponent  macros  in  the  input  file.  (See  ex¬ 
ample  input  file  in  Appendix  C)  The  defcomponent  calls  define  the  basic  geometry  of  the  MCM. 
The  small  heat  sources  found  on  the  chip  surfaces  must  be  defined  as  2nd-subcomponents  by  us¬ 
ing  the  defindsubcomponent  macros.  These  2nd-subcomponents  are  used  to  define  point  heat 
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sources  for  the  first  submodel,  and  the  2nd-subcomponents  in  the  critical  region  are  used  in  the 
second  submodel  analysis.  If  no  2nd-subcomponents  are  defined  and  the  submodel  option  of  run- 
imcma  is  specified,  an  error  will  be  signaled  at  some  point  in  the  analysis. 
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Appendix  C.  Example  of  a  Thermal  Submodeling  Input  File 

Mode;COMMON-LISP;  Package Base:  10  -* 

*-*  File:  /user2/sheehy/imcma2/verif.lisp  *-* 

*-*  Edited-By:  sheehy  *-* 

Last-Edit:  Wed  Apr  13  16:07:41  1994 
*-*  Machine:  hinault.ecs.umass.edu  *-* 


5  5  5  5 

5  5  5  5 

******************************************************************* 

* 

5  5  5  5 

;;;;  *  THERMAL  SUBMODELING  VERIFICATION  EXAMPLE 

....  * 

55  5  5 

5  5  5  5 

4:5f:3|«5j:5H*************************************************  +  ************ 

5  5  55 

:t!SK***************************************************************** 

555 

;;;  Written  by:  Michael  Sheehy 

;;;  UMass  Amherst 

555 

...  :*;************************************* 

5  5  5 

555 

;;;  3-18-94  File  created. 

555 

...  *5j:************************************ 

5  5  5 


(in-package  "IMCMA") 


;;  This  must  come  first: 

(defme-device  "SUBMOD  MCM" 

:filename  "submod-ex" 

:size  (15.00  8.00  1.00)) 

;;  Materials  come  next: 

(defmaterial  :silicon  :tk  .0539  :alpha  0.57e-5) 
(defmaterial  :aluminum  :tk  .15  :alpha  0.27e-5) 


;;  The  substrate: 
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(defcomponent  ;SUBSTRATE-1  rsubstrate 
:size  (15.00  8.00  1.00) 

:  prescribed-temperature-surfaces  ((:bottom  20)) 
:  material  :  aluminum) 

;;  The  chips; 

(defcomponent  :CHIP-1  :chip 
:size  (4.0000  4.000  .25) 

:x  (+  4.00) 

:y  (+  4.00) 

:z(+  1.00) 

:xy-alignment  ;centered 
:power-dissipation  10.341 
imaterial  :silicon) 

(defcomponent  :CHIP-2  ;chip 
:size  (4.0000  4.0000  .25) 

:x  (+11.00) 

:y  (+  4.00) 

;z(+  1.00) 

:xy-alignment  ;centered 
:power-dissipation  5.2 
:  material  :  silicon) 

;;  The  heat  sources; 

(def2ndsubcomponent  ;HS-1  ;heat-source 
;size  (0.5  0.5  0.05) 

;x(+2.75) 

;y  (+2.75) 

;xy-alignment  ;centered 
;power-dissipation  1.5 
;  material  ;  silicon) 

(defZndsubcomponent  ;HS-2  ;heat-source 
;size  (0.5  0.5  0.05) 

;x  (+  4.00) 

;y  (+2.75) 

;xy-alignment  ;centered 
;power-dissipation  1.2 
;material  ; silicon) 

(def2ndsubcomponent  ;HS-3  ;heat-source 
;size  (0.5  0.5  0.05) 
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:x(+5.25) 

:y(+2.75) 

:xy-alignment  icentered 
:power-dissipation  0.7 
imaterial  :silicon) 

(def2ndsubcomponent  :HS-4  :heat-source 
:size  (0.5  0.5  0.05) 

:x(+2.75) 

:y  (+  4.00) 

;xy-alignment  xentered 
:power-dissipation  2.0 
:material  :silicon) 

(def2ndsubcomponent  :HS-5  :heat-source 
:size  (0.5  0.5  0.05) 

:x  (+  4.00) 

:y  (+  4.00) 

:xy-alignment  xentered 
:power-dissipation  1.7 
imaterial  :silicon) 

(def2ndsubcomponent  :HS-6  :heat-source 
:size  (0.5  0.5  0.05) 

:x  (+  5.25) 

;y  (+  4.00) 

:xy-alignment  xentered 
:power-dissipation  0.4 
:material  :silicon) 

(def2ndsubcomponent  :HS-7  :heat-source 
:size  (0.5  0.5  0.05) 

:x  (+  2.75) 

:y  (+  5.25) 

:xy-alignnient  xentered 
:power-dissipation  0.8 
:material  :silicon) 

(def2ndsubcomponent  :HS-8  :heat-source 
:size  (0.5  0.5  0.05) 

;x  (+  4.00) 

:y  (+5.25) 

:xy-alignment  xentered 
:power-dissipation  1.3 
material  :silicon) 
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(defZndsubcomponent  :HS-9  :heat-source 
:size  (0.5  0.5  0.05) 

:x(+5.25) 

;y(+5.25) 

:xy-alignment  xentered 
:power-dissipation  0.741 
imaterial  isilicon) 


;;;  End  of  File 
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Appendix  D.  Implementation  of  Thermal  Submodeling 

The  following  section  details  the  new  files  and  knowledge  sources  created  for  thermal  sub¬ 
modeling,  as  well  as  the  changes  in  the  existing  knowledge  sources.  All  changes  are  documented 
in  the  code,  but  these  descriptions  offer  an  overview  of  the  changes. 

/mcOTfl-system.lisp 

This  file  was  modified  to  include  two  submodeling  knowledge  sources,  model-submodel- 
transition-ks  and  submodel-2nd-submodel-transition-ks,  and  one  additional  file,  imcma- 
interpolation,  in  the  IMCMA  package.  This  ensures  that  the  new  files  are  loaded  with  the  rest  of 
IMCMA. 

imcma-main.lisp 

The  function  run-imcma,  which  starts  the  operation  of  IMCMA,  is  modified  to  contain  ad¬ 
ditional  options  which  are  needed  for  submodeling.  These  options  are  discussed  in  Appendix  B. 
The  init-imcma-internal  function  has  been  modified  to  initialize  some  new  global  variables  cre¬ 
ated  for  the  submodeling  process. 

/mc/Mfl-model-units.lisp 

The  existing  unit  class  definitions  in  this  file  have  been  changed  through  the  addition  of 
new  slots,  and  the  class  names  have  been  modified  to  accommodate  the  storage  of  model  and  sub¬ 
model  unit  classes  on  different  blackboard  spaces.  Additional  submodel  unit  classes  have  also 
been  added  to  this  file. 

/wcma-unit-mappings.lisp 

New  unit  mappings  for  submodel  unit  classes  and  submodel  event  definitions  are  added  to 
this  file. 
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/mcma-preamble.lisp 

Many  global  variables  used  in  submodeling  are  defined  in  imcma-preamble,  and  all  sub¬ 
model  blackboard  and  space  definitions  are  included  in  this  file. 

exodus-utilities.lisp 

The  function  read-analysis-results  has  been  modified  to  accommodate  the  submodel  node 
and  element  unit  classes. 

//Mcma-graphics.lisp 

Display  functions  have  been  modified  to  display  submodel  unit  instances,  and  methods  have 
been  altered  to  accommodate  model  and  submodel  unit  classes,  through  the  specification  of  their 
superclasses. 

/mcma-interpoiation.lisp 

Imcma-interpolation  is  a  new  file,  and  contains  an  interpolation  subroutine  that  is  called  be¬ 
tween  the  levels  of  submodeling.  Also  included  are  functions  needed  for  finding  the  elements  to 
interpolate  from. 


KNOWLEDGE  SOURCES:  The  following  files  contain  IMCMA's  knowledge  sources. 
All  of  the  files  with  the  exception  of  input-model-ks.lisp  have  functions  that  have  been  modified 
to  accommodate  submodel  unit  classes.  Some  functions  have  new  arguments  which  specify  the 
unit  classes  used  at  the  current  level  of  modeling.  These  arguments  are  necessary  to  carry  out  op¬ 
erations  on  the  correct  unit  classes.  In  addition,  global  variables  containing  unit  class  specifica¬ 
tions  are  used  in  retrieval  functions  in  most  or  all  of  the  knowledge  sources.  All  the  files,  again 
with  the  exception  of  input-model-ks,  contain  multiple  knowledge  source  definitions  for  the 
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multiple  levels  of  modeling.  These  KS  definitions  can  be  consolidated  into  one  for  each  file  if 


necessary. 

input-model-ks.lisp 

The  major  modification  to  this  file  is  the  addition  of  the  defsubcomponent  and 
defindsubcomponent  functions,  which  are  currently  called  from  the  IMCMA  input  file  to  define 
unit  instances  of  submodel  components.  Another  small  modifications  is  made  to  input-model-ks 
to  trigger  the  model-submodel-transition-ks,  which  will  be  activated  after  all  other  modeling 
knowledge  sources  are  carried  out. 

complete-model-ks.Iisp 

A  function  has  been  added  to  complete-model-ks  called  change-material-numbers .  This 
function  prevents  exodus  component  block  errors  caused  by  materials  that  are  not  used  at  the  cur¬ 
rent  level  of  modeling.  Another  small  modification  has  been  made  to  trigger  the  correct  knowl¬ 
edge  sources  for  the  current  level  of  modeling. 

generate-mm-regions-ks,iisp 

The  only  addition  is  a  conditional  statement  that  triggers  the  correct  knowledge  source  for 
the  current  level  of  modeling. 

generate-2d-mesh-ks.lisp 

No  modifications  besides  the  changes  described  above  for  all  knowledge  sources. 

extrude-component-ks.lisp 

The  find-wells  function  has  been  modified  to  add  extrusion  layers  for  submodels.  Without 
this  modification,  the  submodels  are  created  with  only  one  layer  of  elements  through  the  thick¬ 
ness.  The  code  can  easily  be  changed  to  obtain  the  desired  number  of  layers. 
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combine-3d-meshes-ks.lisp 

In  the  first  level  of  submodeling,  point  heat  sources  are  created  to  approximate  the  heat 
added  by  the  small  2nd  subcomponent  heat  producing  regions.  The  function  link-prescribed- 
surfaces-to-  elements  has  been  modified  to  identify  the  heat  producing  regions  in  the  submodeling 
region,  and  to  create  point  heat  sources  to  approximate  the  effect  of  these  regions. 

analyze-3d-mesh-ks.lisp 

Some  extensive  additions  have  been  made  to  the  function  invoke-feecap  in  order  to  carry 
out  the  interpolations  necessary  for  prescribing  boundary  conditions  to  a  submodel.  Invoke- 
feecap  has  also  been  changed  to  identify  the  materials  used  at  the  current  level  of  submodel,  and 
to  write  only  the  necessary  materials  to  the  FEECAP  input  file. 

model-submodel-transition-ks.Iisp 

This  is  a  completely  new  knowledge  source  that  performs  many  operations  necessary  for  the 
commencement  of  the  first  level  of  submodeling.  The  display  dimensions  are  changed  for  graphi¬ 
cal  display  of  the  submodel  and  graphics  commands  are  called  to  produce  the  display.  The  region 
to  be  submodeled  is  identified,  and  the  macroscopic  model  elements  which  will  be  used  in  inter¬ 
polation  are  identified  and  stored  in  a  global  variable.  The  submodel  component  unit  instance  is 
created,  and  if  the  adhesive  is  to  be  included  in  the  submodel,  the  die  attach  unit  instance  is  also 
created.  In  the  function  init-imcma-submodel-internal,  several  global  variables  are  set  for  sub¬ 
model  use. 

submodel-2nd-submodel-transition-ks.lisp 

This  knowledge  source  performs  many  of  the  same  functions  as  model-submodel-transition- 
ks,  except  for  the  initiation  of  the  second  submodeling  process.  The  major  difference  is  that  a  die 
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attach  is  not  included  in  the  second  submodel.  Also,  any  2nd  subcomponents  that  exist  on  the 
blackboard  and  are  not  used  in  the  second  submodel  are  identified  and  removed  from  the 
database. 
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Appendix  E.  A  Generic  3-D  Element  Formulation 

The  related  equations  and  a  brief  description  of  the  eight  noded  isoparametric  brick  element 
formulation  for  elastostatic  analysis  is  discussed  in  this  appendix. 

The  Galerkin's  weighted  residual  statement  for  element  e  is  written  as 

j^.N‘{mdr+!^.N‘{R‘a}dn  =  0,  i=l,2,...W'  (47) 

where  is  the  element  surface  area,  is  the  element  volume,  is  the  element  shape  func¬ 
tion  associated  for  the  ith  node  of  the  element,  N‘^  is  the  number  of  nodes  per  element,  and 
{Rr},{Ra}  are  the  boundary  and  domain  residual  vector  functions  respectively  given  by 

{%}=V.[cj‘]  +  {f}  (48) 

{rr}  =  {'">?} -{A},  [o'] 

In  the  above  equations  {«}  is  the  unit  normal  vector,  is  the  component  of  the  externally 

applied  traction  vector  in  the  normal  direction  on  the  element  surface,  {/® }  is  the  body  force  act¬ 
ing  on  the  element  and  [o^]  is  the  element  stress  tensor. 

The  domain  and  boundary  residual  functions  given  by  Equations  (48)  and  (49)  represent  a 
lack  of  satisfaction  of  differential  equilibrium  within  the  element  and  on  the  element  surface.  If 
these  residual  functions  are  identically  equal  to  zero  through  out  the  element  volume  and  every¬ 
where  on  the  element  surface,  then  the  differential  equilibrium  equations  are  completely  satisfied 
in  this  part  of  the  system  domain  and  the  stress  tensor  [a®]  is  the  exact  stress  tensor  admitted  by 
the  theory  of  elasticity. 

Substitution  in  the  Galerkin  weighted  residual  statement  gives 
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‘j’ ({ }-{«}•  [a"])flr  +  jA^^(V«  [o'] +  {/})JQ  =  0,  /=  1,2,...,AA^  (50) 

r‘- 

Integrating  by  parts  the  first  term  of  the  domain  integral  yields 

j Af‘{<”'T‘-}ar-  J  vrr, . [airffi+  / N‘{f}dn  (si) 

r‘'  Q<-  Q.’ 

The  first  term  represents  the  element  nodal  loads  due  to  an  externally  applied  traction  vector 
along  the  element  edges.  Upon  assembly,  this  vector  will  contain  zero  load  values  for  all 
nodes  for  which  there  are  no  externally  applied  traction  vectors  on  interelement  boundaries  and  no 

externally  applied  point  nodal  loads.  The  last  term  represents  the  element  nodal  loads  to  the  body 
force  vector  (f }  .  Body  forces  in  elasticity  are  most  commonly  due  to  gravity  or  inertia  forces  in 
which  a  dynamic  analysis  has  been  converted  into  a  static  analysis  by  D'Alembert's  principle. 

By  letting  i  increment  over  all  nodes  of  the  element  in  Eqn.  (49)  and  using  nontensor  notation, 

we  can  write  the  set  of  element  equations  in  the  form 

c>^ 

a,. 

/  L/Q=K}  +[r}] 

U-- 

where  ^  - 

0  o' 

0  Nly  0 
0  0 
N‘'  0 

i,y  ■'  /,.r  ^ 

0  Ky 
Ny  0  AL  _ 

{r?}  =  ^[A^']"{T}tr 
r*' 

[r}]  =  j  [N^V{f]dn 

Q.‘ 

~  N\  0  0  A^  0  0  ...  Nl,.  0  0 

[A^]  =  0  At  0  0  0  ...  0  N%  0 

0  0  At  0  0  At  ...  0  0  Nl, 


92 


The  term  on  the  left  hand  side  of  the  above  equation  gives  rise  to  the  stiffness  matrix  of  the  ele¬ 


ment.  We  now  introduce  the  constitutive  relations 

{a}  =  [£]({s}-{8}o)  (57) 

The  matrix  [E]  is  and  the  strain  displacement  relations 

"lOOOOOOOO 
000010000 
000000001 
0  1  0  1  0  0  0  0  0 

0  0  0  0  0  1  0  1  0 

0  0  1  0  0  0  1  0  0 

where  u,v,  and  w  are  the  x,  y,  and  z  components  of  the  displacement  field,  respectively  and 

is  the  gradient  operator  with  respect  to  the  x,y,z  coordinate  system. 

For  assumed  displacement  field  based  finite  elements,  the  element  displacement  field 

{u",v",w®}’''  is  interpolated  from  element  nodal  displacements  {u^i,v®i,w^j}’^  using  the  element  shape 

functions  (interpolating  polynomials)  Nj': 

v 

E  NW, 

;=1 

V‘' 

Sacv' 


and  the  matrix  [N*"]  has  been  previously  defined  in  Eqn.  (56).  We  use  Eq.  (58)  to  substitute  for  the 
displacment  field  into  Eq.  (57)  and  obtain 

{e‘'}  =  [B']K}  (60) 
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and  the  matrix  [B']  in  Eqn.  (60)  is  identical  to  the  [B"]  matrix  defined  in  Eqs.  (53)  and  (54).  Using 
Eqn.  (60Eqn.  (57),  the  following  expression  is  obtained  for  the  element  stress  vector. 

{cr'}  =  [£']([B']{^/'}-{eo})  (61) 

At  this  point  we  can  now  substitute  for  the  stress  vector  given  by  Eqn.  (61)  into  Eqn.  (52)  to 


obtain 


J  -  {eo})  dn  =  {K}  +  {r}} 


Rearranging  this  expression  produces  the  final  element  equations: 


where 


WV,d‘}  =  {r\}  +U)  +  {r».} 

[K‘\  =  j  [B’YWWB’YCI 

{'■«}  = 


All  that  remains  to  do  is  to  adopt  the  isoparametric  formulation,  since  it  makes  the  integrals 
easier  to  evaluate  numerically,  leads  to  compatible  elements,  and  makes  the  element  shape  func¬ 
tions  independent  of  the  element  geometry. 

Since  the  shape  functions  are  written  in  terms  of  the  natural  coordinates  ^ ,  p  and  p  we 
must  convert  the  gradient  operator  with  respect  to  x,  y  and  z  to  the  gradient  operator  with  respect 


to  the  natural  coordinates; 


where 


[ri  =  m- 


=  [r®]^^pv 

•^w  '^n 

-  '^12  ^22 
-^1  *^2  *^3 
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and 


'^W 

= 

54 

Omrt' 

^  5V,-  g 
-  Zj  g. 

/=! 

= 

54 

^  5^*'  g 

dz‘ 

54 

/=l 

A\ 

= 

5n 

y  5^,.  g 

“  ^  dt]^‘ 
;=! 

fj2 

= 

Sy" 

5n 

=  ^ 

i=\ 

6^ 

5n 

"  h 

A\ 

= 

5p 

_y 

""y  5p^< 
/=1 

Ai 

= 

5p 

=  Z— / 

Ai 

= 

dz" 

dp 

_y 
(=1  ^ 

(68) 


[J®]  is  the  3x3  Jacobian  transformation  matrix  from  the  x,y,z  system  to  the  natural  coordinate 
system. 

Most  of  the  formulation  discussed  above  contains  equations  that  are  similar  to  those  of  the 
tetrahedral  element  formulation  (Appendix  1).  The  equations  used  in  this  appendix  are  generic 
equations  and  may  be  used  for  any  3-D  element  type. 
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Appendix  F.  Determining  the  Critical  Interconnect 

/MCM4-interconnect.lisp: 

This  module  was  developed  to  define  the  unit-class  for  interconnects.  The  various  slots  for 
the  unit-class  interconnect  were 

•  type  -  indicating  type  of  interconnect  (e.g.:  wirebond,  TAB  etc.)  material 

•  indicating  type  of  material  of  interconnect. 

•  start-point  •  location  of  start-point  of  interconnect  on  the  MCM. 

•  end-point  -  location  of  end-point  of  interconnect  on  the  MCM. 

•  start-displacements  -  stores  the  displacements  (interpolated  from  the  macroscopic  mesh)  at 
the  start-point  of  the  interconnect 

®  end-displacements  -  stores  the  displacements  (interpolated  from  the  macroscopic  mesh)  at 
the  end-point  of  interconnect. 

•  static-temperatures  -  stores  the  temperatures  (interpolated  from  the  macroscopic  mesh)  at 
the  start-point  of  interconnect  end-ic-temperatures 

•  stores  the  temperatures  (interpolated  from  the  macroscopic  mesh)  at  the  end-point  of 
interconnect 

•  relative-displacement  -  stores  the  relative  displacement  of  the  interconnect. 

input-model-ks.Iisp: 

This  module  was  modified  to  include  the  defmacro  for  the  interconnect  according  to  the 
unit-class  definition.  The  defmacro  "defintercomiect  "  enables  the  user  to  input  the  attributes  of 
each  interconnect  and  correspondingly  each  instance  of  the  interconnect  is  saved  onto  the 
blackboard. 

/MCM4-preamble.lisp  : 
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The  blackboard  path  for  each  instance  of  the  interconnect  is  defined  and  space  is  created  to 
store  the  instances  of  the  interconnect.  The  path  for  the  interconnect  instance  is  "model  3d-model 
interconnects".  Any  information  regarding  the  intereonnects  is  retrieved  by  specifying  this  path. 

feecap.f : 

Feecap  had  to  be  modified  to  write  out  the  nodal  displacements  in  the  exodus  format 

/MCMd-model-units.Iisp  ; 

New  slots  had  to  be  defined  for  the  nodal-displacements  unit  class  to  store  the  displace¬ 
ments  at  each  node.  Moreover,  since  a  combined  analysis  has  to  be  performed  for  the  identifica¬ 
tion  of  critical  interconnect  a  unit  class  POINT-DISPLACEMENT-CONSTRAINT  was  added. 
Thus  while  defining  the  components  sueh  as  chip  and  substrate  in  the  input  file,  point  displace¬ 
ment  constraints  could  be  specified  on  the  macroscopic  model  as  (lower-left- front-comer  (10  0)), 
where  1  indicates  constraint  and  0  indicates  free.  Hence,  enabling  displacement  constraints  to  be 
specified,  combined  analysis  could  thus  be-carried  out. 

exodus-utilities.Iisp  : 

This  module  had  to  be  modified  to  read  the  displacements  of  each  node  written  in  the  exo¬ 
dus  format.  The  displacements  of  each  node  are  then  stored  onto  the  blackboard. 

/MCM4-unit-mappings.lisp 

In  this  module  an  event  is  defined  to  enable  to  trigger  the  critical  interconnect  knowledge 
source. 

critical-interconnect-ks.lisp; 

This  knowledge  source  was  developed  to  identify  the  critical  interconnect.  Once  the  macro¬ 
scopic  mesh  analysis  is  performed  the  nodal  displacements  of  each  element  is  available  and  can 
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be  retrieved  from  the  blackboard.  The  location  of  each  interconnect  is  identified  with  the  top  ele¬ 
ments  of  the  macroscopic  mmdel  and  displacements  at  the  nodes  of  the  element  which  encom¬ 
passes  the  start  point  or  the  end  point  of  the  interconnect  are  interpolated  onto  the  start  or  end 
point  of  the  interconnect  respectively.  The  interpolation  routine  is  discussed  below 

Interpolation:  The  MCM  is  modeled  with  rectilinear  eight-noded  brick  elements.  A  two- 
dimensional  interpolation  is  carried  out  for  the  interpolation  of  displacem.ents  from  the  nodes  of 
the  element  to  the  start  or  end  point  of  the  interconnect  since  the  start  or  the  end  point  of  the  inter¬ 
connect  lie  on  a  face  of  the  encompassing  element  and  not  in  the  interior  of  an  element.  This  sim¬ 
plifies  our  interpolation  from  three  dimensional  to  two  dimensional.  The  shape  functions  or 
interpolating  functions  are  calculated  for  the  four  nodes  of  a  face  of  the  element  in  terms  of  the 
isoparametric  coordinates.  The  shape  functions  are  given  as: 

^2(4,  h)  =  +^)(i  -"n)  (69) 

=  i(l+^)(l+Tl) 

As  the  real  space  coordinates  of  the  interpolating  point  are  known,  these  shape  functions  are 

then  substituted  in  terms  of  the  isoparametric  coordinates  in  the  following  equations: 

X]  —  N\X\  -f  A^22r2  +  A^3X3 -f- A^4X4  (70) 

yi  —  A^iy'i  +  A^2T2  +  ■^4T4 

where  x,  and  y,  are  the  real  space  coordinates  of  the  interpolating  point  /,  N,  to  N4  are  the  shape 

functions  of  the  element  which  encloses  the  interpolating  point,  and  x,  to  X4,  y,  to  y4  are  the  x  and 

y  coordinates  of  the  nodes  of  the  face  of  the  element  respectively.  The  isoparametric  coordinates 
p  are  then  solved  knowing  the  coordinates  of  the  start  point  or  end  point  of  the  interconnect 
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from  the  two  functions 


4x/-(xi  +X2  +X3  +X4)  =  (-Xi  +X2  +X3  -X4)^  +  (-Xi  -X2  +X3  +X4)ri  +  (xi  -X2  +X3  -X4)^ri 

4v/-Oi  +y2  +y3  +3^4)  =  (-3^1  +y2  +y3  -y4)^  +  (-yi  -y2  +3^3  +3^4)11  +  (yi  -3^2  +3^3  -3M)^h 

(71) 

The  root  (^/,ri/)  which  satisfies  Eqn.  (71)  and  the  condition 

-l^^^l  (72) 

-1  <  T)  <  1 

is  substituted  back  into  the  shape  functions  and  the  displacements  at  point  /  (i.e.  either  the  start 

or  end  point  of  the  interconnect)  is  interpolated  from  nodal  displacements  using  the  formula 

u/  =  Ni(^/,ri/)ui+N2(^/,r]/)u2+N3(^/,ri/)u3+N4(^/,ri/)u4 
V/  =  Ni(^/,r]/)vi  +N2(^/,ri/)v2+N3(^/,ri/)v3+N4(^/,r]/)V4 
m  =  N](^/,ri/)wi-i-N2(^/,ri/)w2+N3(^/,r]/)w3+N4(^/,r\/)w4 


Once  the  displacements  are  interpolated  to  the  start  and  end  point  of  each  interconnect  the 
relative  displacement  is  calculated  for  each  interconnect.  The  relative  displacement  between  the 

start  and  end  points  of  the  interconnect  is  determined  by: 

Ad  =  J(Us  -  Ue^  +  (v^  -  Vef  +  (w^  -  Wef  (7  4) 

The  interconnect  with  the  largest  relative  displacement  is  identified  and  stored  on  the  black¬ 
board  while  the  remaining  interconnect  instances  are  deleted  from  the  blackboard.  Further  sub¬ 
modeling  of  the  critical  interconnect  can  be  performed  by  retrieving  the  instance  from  the 
blackboard. 

/MCM4-main.Iisp  : 
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This  module  has  been  modified  to  include  an  option  for  identifying  the  critical  interconnect. 
While  invoking  xun-imcma  with  the  option  -  interconnect  true  -  the 
CRITICAL-INTERCONNECT-KS.LISP  is  triggered. 

With  these  new  additions  and  modifications  the  identification  of  a  critical  interconnect  is 
fully  functional. 
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Appendix  G.  Modeling  Tips  for  the  Interconnect 

The  user  is  free  to  use  the  wide  variety  of  features  available  in  PRO-ENGINEER  for  model¬ 
ing  the  interconnect.  A  few  methods  used  by  the  developers  are  suggested  below: 

Method  1 

•  Create  default  datum  planes  and  coordinate  axes. 

•  Read  in  the  datum  points  from  the  points.dat  file,  which  is  automatically  created  during  the 
execution  of  the  PRO-MODEL- WRITE-KS. 

•  Use  the  solid  pipe  feature  (for  wirebonds)  to  create  the  interconnect  joining  the  points  us¬ 
ing  a  spline  curve. 

•  Create  the  thermocompression  (for  wirebonds)  on  one  end  face  of  the  interconnect  as  a 
protrusion  feature.  Rounding  all  edges  is  a  good  idea. 

•  Copy  this  feature  on  the  other  end  face  by  specifying  alternate  datum  planes  and  points  us¬ 
ing  the  copy  feature  option  in  PRO-E. 

•  Enter  the  FEM  module  of  PRO-E  and  create  the  tetrahedral  mesh  using  the  tet-mesh  option. 
Do  not  use  the  default  global  element  length  for  the  meshing  process. 

•  Output  the  mesh  to  an  ANSYS  file  bond.ans  choosing  the  linear  structural  element  type 
from  the  menu. 

•  Exit  PRO-E. 

Method  2 

•  Repeat  first  two  steps  from  method  1 . 

•  Use  the  protrusion-sweep  menu  to  create  the  interconnect. 
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•  First  sketch  the  spline  curve  starting  with  the  start  point  and  ending  with  the  end  point  to 
define  the  profile  of  the  interconnect.  Dimension  the  spline  with  reference  to  the  default  da¬ 
tum  planes. 

•  Sweep  the  section  of  the  interconnect  along  the  profile. 

•  Repeat  steps  4  through  8  of  method  1 . 

Note:  While  any  method  may  be  used  to  model  the  interconnect,  the  following  is 
mandatory: 

•  The  start  and  end  points  always  lie  at  the  "z-level"  at  which  the  nodes,  which  are  pre¬ 
scribed  displacements,  lie. 

•  The  start  point  and  end  point  must  be  read  from  the  points.dat  file.  The  choice  of  the  inter¬ 
mediate  points  is  left  to  the  user. 

•  The  start  point  is  the  location  on  the  interconnect  where  the  section  of  the  interconnect 
starts.  The  flattened  region/the  region  where  the  interconnect  is  fixed  to  the  chip  or  sub¬ 
strate  must  extend  away  from  this  section.  The  same  rule  applies  to  the  end  point.  The 
length  of  this  region  is  the  thermocompression  length. 
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Appendix  H.  Interconnect  Modeling  and  Analysis  -  Knowledge  Sources 
Added  and  Modified 

The  following  knowledge  sources  have  been  added  to  IMCMA  for  the  interconnect  model¬ 
ing  and  analysis: 

(1)  pro-model-ks.lisp  :  This  knowledge  source  is  responsible  for  interfacing  IMCMA  with 
PRO-ENGINEER.  When  executed,  it  opens  up  a  PRO-ENGINEER  window.  Some  user 
instructions  are  also  output  on  the  LISP  window. 

(2)  tetra-3d-mesh-ks.lisp:  This  knowledge  source  reads  in  the  mesh  data  from  the  ANSYS 
file  and  instantiates  nodes,  elements  on  the  blackboard  and  also  links  the  interconnect  in¬ 
stance  with  the  mesh  data. 

(3)  tetra-mesh-analyze-ks.lisp:  This  knowledge  source  determines  the  nodes  that  lie  at  the 
start  and  end  of  the  interconnect  and  assigns  the  corresponding  displacements  to  these 
nodes.  Further  it  accesses  the  mesh  data,  material  information  and  boundary  conditions 
from  the  blackboard  and  writes  out  the  FEECAP  input  file. 

(4)  read-results-ks.lisp:  This  knowledge  source  reads  in  the  post  analysis  results  obtained 
from  FEECAP  and  the  results  are  written  back  to  the  blackboard.  The  data  read  includes 
nodal  displacements  and  stresses. 

The  following  knowledge  sources  have  been  modified  in  the  IMCMA  system  to  integrate 
the  modeling  of  the  interconnect  with  the  other  existing  modules: 

(1)  imcma-systenulisp:  An  addition  has  been  made  to  the  list  of  knowledge  sources  to  in¬ 
clude  all  the  new  knowledge  sources  written  for  interconnect  modeling. 
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(2)  imcma-main.lisp:  A  counter  for  the  number  of  interconnect  nodes  and  elements  has  been 
introdueed. 

(3)  imcma-preamble.lisp:  This  knowledge  source  is  modified  to  inelude  blackboard  spaces 
for  the  interconnect,  interconnect-3 d-node  and  interconnect-3 d-element  and  also  to  define 
the  global  variable  for  the  number  of  nodes  with  preseribed  displaeements. 

(4)  imcma-unit-mappings.lisp:  The  events  and  signals  for  the  new  knowledge  sources  have 
been  added  to  the  list  of  existing  events. 
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Appendix  I.  Tetrahedral  Element  Formulation  -  A  Summary  of  Formulas 

The  finite  element  code  for  linear  tetrahedral  formulation  has  been  written  using  references 

[10]  -  [12],  Reference  [10]  has  been  utilized  for  the  thermal  analysis  formulation  while  reference 

[11]  has  been  used  to  formulate  the  static  analysis  code.  The  general  formula  for  the  element  ther¬ 
mal  conductance  matrix  with  coordinate  axis  aligned  with  orthotropic  material  conductivity  tensor 


,  r  f  dN]  bn; 


+  K‘yy 


dy  dy  oz  oz 


ij=  1,2, 


where  Q"  is  the  element  volume,  K^yy,  and  are  the  element  thermal  conductivities  in  the 
x,y,  and  z  directions,  respectively,  N‘ 's  are  the  element  nodal  shape  functions,  and  N^’  is  the  num¬ 
ber  of  nodes  for  element  e.  For  the  four-noded  tetrahedral  element,  Eqn.  (75)  can  be  integrated  in 
closed  form  to  give 

b\b\  b\b2  bxbz  b\bi^  ciCi  C]C2  C1C3  C1C4 

,  ,  b2b\  b2b2  ^2^3  ^2^4  K.'yy  C2C\  C2C2  C2C3  C2C4 

^  ^  36Q^  b^bx  b3b2  b^b^  63^4  360.^  c^cx  C3C2  C3C3  C3C4 

b^bx  ^4^2  ^4^3  b/^b/\  C4C1  C4C2  C4C3  C4C4 

dxdx  dxd2  dxd^  did 4 
Kzz  d2dx  d2d2  d2d3  J2<^4  (76) 

36^2^  d^dx  d3d2  d^d^  d^d^ 
d^dx  d4d2  dixd^  d/[d4 

where 

bi  =  (y2-  y4)(z3-  Z4)  -  (y3  -  y4)(z2-  Z4) 
b2  =  (y3-  y4)(zr  Z4)  -  (y,  -  y4)(z3-  Z4) 

b3  =  (yr  y4)(z2-  Z4)  -  (y2  -  y4)(zi-  Z4) 

b4=  -  (b|  +  b2  +  b3) 


105 


C,  =  (X3-  X4)(Z2-  Z4)  -  (X2  -  X4)(Z3-  Z4) 

C2  =  (X,-  X4)(Z3-  Z4)  -  (X3  -  X4)(Z,-  Z4) 

C3  =  (X2-  X4)(Z,-  Z4)  -  (X,  -  X4)(Z2-  Z4)  (77) 

C4=  -(C,+  C2  +  C3) 

d,  =  (X2-  X4)(y3-  y,)  -  (X3  -  X4)(y2-  y4) 

d2  =  (X3-  X4)(y,-  y4)  -  (x,  -  X4)(y3-  y4) 
d,  =  (x,-  X4)(y2-  y4)  -  (X2  -  X4)(y|-  y4) 
d4  =  -(d,  +d2  +  d3) 

The  elemental  K  matrix  for  a  thermal  analysis  for  a  four  noded  tetrahedral  element  is  a  4  X 
4  matrix  (N^  ■  NDF  x  ■  NDF)  where  NDF  is  the  number  of  degrees  of  freedom  per  node 
whieh  is  unity  for  thermal  problems.  In  case  of  a  static  analysis  the  elemental  K  matrix  is  a  12  X 

12  matrix  since  there  are  three  degrees  of  freedom  per  node  corresponding  to  displacements  in  the 
x,y,  and  z  direetion. 

For  a  thermal  analysis  the  elemental  load  vector  is  as  follows 

{'■'}  =  {4>  +  {''e}  +  (4}  (78) 

where  {r\^}  is  the  thermal  load  due  to  prescribed  face  flux  q^,  {Fq}  is  the  load  vector  due  to  in¬ 


ternal  heat  generation,  and  {/‘^jis  the  thermal  load  due  to  convection  boundary  conditions  on  a 
face.  The  integral  expressions  and  resulting  values  for  these  vectors  are 


123 

1110 

1111 

4  L 

_  FFuA 
~  3 

^[011 

where  we  have  taken  Q^,h  ,  and  T%  to  be  constant  over  the  element  integrals.  A, 23  and  A2U 
are  the  surface  areas  of  element  faces  defined  by  nodes  1-2-3  and  2-3-4,  respectively,  and  we 

have  assumed  the  flux  to  be  applied  on  face  1-2-3  and  convection  boundary  conditions  acting  on 
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face  2-3-4  for  illustrative  purposes.  The  above  formulae  show  that  all  integrals  are  reduced  to 
straightforward  calculations  and  evaluations  at  nodal  points. 

In  case  of  a  static  analysis,  the  strain  vector  for  an  element  is  defined  as 


{£«}  =  [ 


Sx  Cz 


Jxy  Jyz  yyz  } 


or  in  terms  of  the  [B]  matrix 


where 


bi  0  0 

0  Ci  0 

r^e] - L_,  ^  ^  >  /  =  1  2  3  4 

0  di  Ci 
di  0  bi 

< 

The  b,  c  and  d's  are  same  as  previously  defined  and  Q*  is  the  volume  of  the  element. 
In  general  the  stress  vector  in  three-dimensional  elasticity  is  given  by 


=  [£‘’]({s«)-{ej})  +  {a'} 


where  [E*^]  is  the  constitutive  material  matrix,  {8*^0}  and  {a"o}  are  the  initial  strain  and  stress 
acting  on  the  element.  Initial  stresses  may  be  induced  by  manufacturing  processes  as  residual 
stresses  in  the  material  and  are  difficult  to  measure  or  predict.  Initial  strains  are  typically  due  to 
thermal  expansion  or  contraction  of  the  material  due  to  heating  or  cooling  of  the  material  from  its 
strain  free  reference  temperature  Thermally  induced  initial  strains  are  hydrostatic  in  nature 
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and  are  given  by 


0 
0 
0 

where  T  is  the  element's  temperature  field  and  a%^,  and  are  the  elements  coefficients  of 
thermal  expansion  in  the  x,y,  and  z  direction,  respectively.  The  coupling  of  thermal  behavior  with 
stress  analysis  is  now  evident,  since  the  temperature  field  within  an  element  is  required  to  com¬ 
pute  initial  strains.  The  initial  strain  vector  produces  a  static  load  vector  according  to  Eqn.  (65) 
and  affects  the  stress  state  according  to  Eqn.  (85). 
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Appendix  J.  Sample  Input  File  for  FEECAP  2.5. 


The  following  is  the  INPUT.DAT  file  for  the  new  version  of  FEECAP.  As  ean  be  noted 
from  the  analysis  type  index  (ANAL  TYP  =  3),  the  input  file  is  for  a  combined  analysis.  The  lines 
marked  '@1'  are  absent  for  a  static  analysis  and  those  marked  '@2'  are  absent  for  a  thermal 
analysis. 

The  model  is  that  of  a  cube  the  vertices  of  which  are  nodes.  There  is  an  additional  node  at 


the  center  of  the  cube.  The  model  has  nine  nodes.  Twelve,  four  noded  tetrahedral  elements  are  re¬ 


quired  to  mesh  this  model.  The  elemental  data  in  the  input  file  shows  the  nodal  connectivity  for 

each  element.  All  possible  boundary  conditions  that  may  exist  have  been  shown. 

TITLE 

TETRAHEDRAL  ELEMENT  INCORPORATION 
NUMNP,  NUMEG,  MODEX,  ANAL.  TYP,  IDTYP 
9,  1,  1,  3,  3 


NODAL  DATA 


1 

1 

1 

1 

0 

0.0000000 

1.000000 

1.000000 

2.000000 

2 

0 

1 

1 

0 

0.0000000 

-1.000000 

1.000000 

2.000000 

3 

0 

0 

1 

0 

0.0000000 

-1.000000 

-1.000000 

2.000000 

4 

0 

0 

1 

1 

40.000000 

1.000000 

-1.000000 

2.000000 

5 

0 

0 

1 

0 

0.0000000 

1.000000 

1.000000 

0.000000 

6 

0 

0 

0 

1 

80.000000 

-1.000000 

1.000000 

0.000000 

7 

0 

1 

0 

0 

0.0000000 

-1.000000 

-1.000000 

0.000000 

8 

0 

0 

0 

1 

90.000000 

1.000000 

-1.000000 

0.000000 

9 

0 

0 

0 

0 

0.000000 

0.000000 

0.000000 

1.000000 

NELTYP,  #  ELEM,  #  MAT'LS,  NPE 
3,  12,  1,  4 

MAT'LNO.,  E,  NU,  MIN_TOL,  MAX_TOL,  BLOCK# 

1,  0.540E+08,  0.2200,  0.1000,  0.2000, 1 
TKr,  TKs  ,  TKz,  BETA  [DEGREES] 

10.00000,  10.00000,  10.00000,  0.000000 

Alpha-R,  Alpha-S,  Alpha-Z  Tref 

0.640000007E-05,  0.640000007E-05,  0.640000007E-05,  0.0000000000 
NODAL  CONNECTIVITY 


@1 

(^1 

@1 

@1 
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1  1  5  4  9  1  0 

2  8  4  5  9  1  0 

3  7  4  8  9  1  100 

4  7  3  4  9  1  0 

5  2  4  3  9  1  0 

6  2  1  4  9  1  0 

I  2  3  1  9  \  100 

8  2  7  6  9  1  0 

9  2  5  1  9  1  0 

10  2  6  5  9  1  0 

II  7  8  6  9  1  100 

12  6  8  5  9  1  0 


FACE  CONVECTION  CARDS 
2 

ELEM#  FACE#  H  TAMB 
2  1  0.6  100 
4  4  0.8  100 


@1 

@1 

@1 

@1 

@1 


FACE  FLUX  CARDS 
2 

ELEM#  FACE#  FLUX 
6  1  100 
3  2  100 


@1 

@1 

@1 

@1 

@1 


THERMAL  LOAD  DATA 
2 

2  100 
4  100 


@1 

@1 

@1 

@1 


STATIC  LOAD  DATA 

4 

5  3  -100 

6  3  -100 

7  3  -100 

8  3  -100 


@2 

@2 

@2 

@2 

@2 

@2 


PRES  NODAL  DISP 
2 

5  3  le-5 
7  2  2e-5 


*U.S.  GOVERNMENT  PRINTING  OFFICE:  1995-610-126-50129 
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MISSION 

OF 

ROME  LABORATORY 


Mission.  The  mission  of  Rome  Laboratory  is  to  advance  the  science  and 
technologies  of  command,  control,  communications  and  intelligence  and  to 
transition  them  into  systems  to  meet  customer  needs.  To  achieve  this, 
Rome  Lab: 


a.  Conducts  vigorous  research,  development  and  test  programs  in  all 
applicable  technologies; 

b.  Transitions  technology  to  current  and  future  systems  to  improve 
operational  capability,  readiness,  and  supportability; 

c.  Provides  a  full  range  of  technical  support  to  Air  Force  Materiel 
Command  product  centers  and  other  Air  Force  organizations: 

d.  Promotes  transfer  of  technology  to  the  private  sector; 

e.  Maintains  leading  edge  technological  expertise  in  the  areas  of 
surveillance,  communications,  command  and  control,  intelligence,  reliability 
science,  electro-magnetic  technology,  photonics,  signal  processing,  and 
computational  science. 


The  thrust  areas  of  technical  competence  include:  Surveillance, 
Communications,  Command  and  Control,  Intelligence,  Signal  Processing, 
Computer  Science  and  Technology,  Electromagnetic  Technology, 
Photonics  and  Reliability  Sciences. 


